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Summary. The first dwarf galaxies, which constitute the building blocks of the 
collapsed objects we find today in the Universe, had formed hundreds of millions of 
years after the big bang. This pedagogical review describes the early growth of their 
small-amplitude seed fluctuations from the epoch of inflation through dark matter 
decoupling and matter-radiation equality, to the final collapse and fragmentation of 
the dark matter on all mass scales above ~ 1O~*M0. The condensation of baryons 
into halos in the mass range of ~ lO^-lO^^M© led to the formation of the first stars 
and the re-ionization of the cold hydrogen gas, left over from the big bang. The 
production of heavy elements by the first stars started the metal enrichment process 
that eventually led to the formation of rocky planets and life. 

A wide variety of instruments currently under design [including large-aperture 
infrared telescopes on the ground or in space {JWST), and low-frequency arrays for 
the detection of redshifted 21cm radiation], will establish better understanding of the 
first sources of light during an epoch in cosmic history that was largely unexplored so 
far. Numerical simulations of reionization are computationally challenging, as they 
require radiative transfer across large cosmological volumes as well as sufficently high 
resolution to identify the sources of the ionizing radiation. The technological chal- 
lenges for observations and the computational challenges for numerical simulations, 
will motivate intense work in this field over the coming decade. 
Disclaimer: This review was written as an introductory text for a series of lec- 
tures at the SAAS-FEE 2006 winter school, and so it includes a limited sample of 
references on each subject. It does not intend to provide a comprehensive list of all 
up-to-date references on the topics under discussion, hut rather to raise the interest 
of beginning graduate students in the related literature. 



1 Opening Remarks 

When I open the daily newspaper as part of my morning routine, I often 
see lengthy descriptions of conflicts between people on borders, properties, 
or liberties. Today's news is often forgotten a few days later. But when one 
opens ancient texts that have appealed to a broad audience over a longer 
period of time, such as the Bible, what does one often find in the opening 
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chapter?... a discussion of how the constituents of the Universe (including 

light, stars and life) were created. Although humans arc often occupied with 
mundane problems, they are curious about the big picture. As citizens of the 
Universe, we cannot help but wonder how the first sources of light formed, 
how life came to existence, and whether we are alone as intelligent beings in 
this vast space. As astronomers in the twenty first century, we are uniquely 
positioned to answer these big questions with scientific instruments and a 
quantitative methodology. In this pedagogical review, intended for students 
preparing to specialize in cosmology, I will describe current ideas about one 
of these topics: the appearance of the first sources of light and their influence 
on the surrounding Universe. This topic is one of the most active frontiers 
in present-day cosmology. As such it is an excellent area for a PhD thesis 
of a graduate student interested in cosmology. I will therefore highlight the 
unsolved questions in this field as much as the bits we understand. 

2 Excavating the Universe for Clues About Its History 

When we look at our image reflected off a mirror at a distance of 1 meter, 
we see the way we looked 6 nano-seconds ago, the light travel time to the 
mirror and back. If the mirror is spaced 10^^ cm = 3pc away, we will see 
the way we looked twenty one years ago. Light propagates at a finite speed, 
and so by observing distant regions, we are able to see how the Universe 
looked like in the past, a light travel time ago. The statistical homogeneity 
of the Universe on large scales guarantees that what we see far away is a fair 
statistical representation of the conditions that were present in in our region 
of the Universe a long time ago. 

This fortunate situation makes cosmology an empirical science. We do not 
need to guess how the Universe evolved. Using telescopes we can simply see 
the way it appeared at earlier cosmic times. Since a greater distance means a 
fainter flux from a source of a flxed luminosity, the observation of the earliest 
sources of light requires the development of sensitive instruments and poses 
challenges to observers. 

We can in principle image the Universe only if it is transparent. Earlier 
than 0.4 million years after the big bang, the cosmic plasma was ionized and 
the Universe was opaque to Thomson scattering by the dense gas of free 
electrons that filled it. Thus, telescopes cannot be used to image the infant 
Universe at earlier times (or redshifts > 10'^). The earliest possible image of 
the Universe was recorded by COBE and WMAP (see Fig. 2). 
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Cosmic-Archeology 



Early time 




can trace the history of the universe. 



Fig. 1. Cosmology is like archeology. The deeper one looks, the older is the layer 
that one is revealing, owing to the finite propagation speed of light. 




Fig. 2. Images of the Universe shortly after it became transparent, taken by 
the COBE and WMAP satellites (see http://map.gsfc.nasa.gov/ for details). The 
slight density inhomogeneties in the otherwise uniform Universe, imprinted hot 
and cold brightness map of the cosmic microwave background. The existence 
of these anisotropies was predicted three decades before the technology for tak- 
ing this image became available in a number of theoretical papers, including 
[355, 308, 297, 338, 282]. 
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Fig. 3. The optical depth of the Universe to electron scattering (upper panel) and 
the ionization fraction (lower panel) as a function of redshift before reionization. Ob- 
servatories of electromagnetic radiation cannot image the opaque Universe beyond 
a redshift of 2: ~ 1100. 

3 Bakground Cosmological Model 
3.1 The Expeinding Universe 

The modern physical description of the Universe as a whole can be traced 
back to Einstein, who argued theoretically for the so-caUed "cosmological 
principle" : that the distribution of matter and energy must be homogeneous 
and isotropic on the largest scales. Today isotropy is well established (see the 
review by Wu, Lahav, & Rees 1999 [389]) for the distribution of faint radio 
sources, optically-selected galaxies, the X-ray background, and most impor- 
tantly the cosmic microwave background (hereafter, CMB; see, e.g., Bennett et 
al. 1996 [36]). The constraints on homogeneity are less strict, but a cosmolog- 
ical model in which the Universe is isotropic but significantly inhomogeneous 
in spherical shells around our special location, is also excluded [155]. 

In General Relativity, the metric for a space which is spatially homoge- 
neous and isotropic is the Friedman- Robertson- Walker metric, which can be 
written in the form 



where a(t) is the cosmic scale factor which describes expansion in time, and 
(i?, Q, (j)) are spherical comoving coordinates. The constant k determines the 
geometry of the metric; it is positive in a closed Universe, zero in a flat Uni- 
verse, and negative in an open Universe. Observers at rest remain at rest, at 
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fixed {R, 0, (j)), with their physical separation increasing with time in propor- 
tion to a(t). A given observer sees a nearby observer at physical distance D 
receding at the Hubble velocity H{t)D, where the Hubble constant at time t 
is H{t) = da{t)/dt. Light emitted by a source at time t is observed at f = 
with a redshift z = l/a(i) — 1, where we set a(t = 0) = 1 for convenience (but 
note that old textbooks may use a different convention). 

The Einstein field equations of General Relativity yield the Priedmann 
equation (e.g., Weinberg 1972 [376]; Kolb & Turner 1990 [205]) 

i^^W = — P-^, (2) 

which relates the expansion of the Universe to its matter-energy content. For 
each component of the energy density p, with an equation of state p = p{p), 
the density p varies with a{t) according to the equation of energy conservation 

d{pR^) = -pd{R^) . (3) 

With the critical density 

defined as the density needed for A; = 0, we define the ratio of the total density 
to the critical density as 

= . (5) 

PC 

With fim, a-nd fir denoting the present contributions to fi from matter 
(including cold dark matter as well as a contribution Q}, from baryons), vac- 
uum density (cosmological constant), and radiation, respectively, the Fried- 
mann equation becomes 



Hit) 



(6) 



where we define Hq and = f2m, + f^A + to be the present values of H 
and f2, respectively, and we let 

k 
Hq 

In the particularly simple Einstein-de Sitter model (/?„ = 1, l7yi = = 

fik = 0), the scale factor varies as a{t) cx i^/^. Even models with non-zero 
Qa or approach the Einstein-de Sitter behavior at high redshift, i.e. when 
[1 + z)~:^ \Qm~^ — Ij (as long as Qr can be neglected). In this high-^ regime 
the age of the Universe is, 

(1 + ^)-'/' . (8) 
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The Friedmann equation implies that models with i?^ = converge to the 
Einstein-dc Sitter limit faster than do open models. 

In the standard hot Big Bang model, the Universe is initially hot and the 
energy density is dominated by radiation. The transition to matter domina- 
tion occurs at z ^ 3500, but the Universe remains hot enough that the gas 
is ionized, and electron-photon scattering effectively couples the matter and 
radiation. At 2; ~ 1100 the temperature drops below ~ 3000K and protons 
and electrons recombine to form neutral hydrogen. The photons then decouple 
and travel freely until the present, when they are observed as the CMB [348]. 

3.2 Composition of the Universe 

According to the standard cosmological model, the Universe started at the big 
bang about 14 billion years ago. During an early epoch of accelerated superlu- 
minal expansion, called inflation, a region of microscopic size was stretched to 
a scale much bigger than the visible Universe and our local geometry became 
flat. At the same time, primordial density fluctuations were generated out 
of quantum mechanical fluctuations of the vacuum. These inhomogeneities 
seeded the formation of present-day structure through the process of gravita- 
tional instability. The mass density of ordinary (baryonic) matter makes up 
only a flfth of the matter that led to the emergence of structure and the rest 
is the form of an unknown dark matter component. Recently, the Universe 
entered a new phase of accelerated expansion due to the dominance of some 
dark vacuum energy density over the ever rarefying matter density. 
The basic question that cosmology attempts to answer is: 
What are the ingredients (composition and initial conditions) of the 
Universe cind what processes generated the observed structures in 
it? 

In detail, we would like to know: 

(a) Did inflation occur and when? If so, what drove it and how did it end? 

(b) What is the nature of of the dark energy and how does it change over time 
and space? 

(c) What is the nature of the dark matter and how did it regulate the evolution 
of structure in the Universe? 

Before hydrogen recombined, the Universe was opaque to electromagnetic 
radiation, precluding any possibility for direct imaging of its evolution. The 
only way to probe inflation is through the fossil record that it left behind 
in the form of density perturbations and gravitational waves. Following in- 
flation, the Universe went through several other milestones which left a de- 
tectable record. These include: baryogenesis (which resulted in the observed 
asymmetry between matter and anti-matter), the electroweak phase transition 
(during which the symmetry between electromagnetic and weak interactions 
was broken), the QCD phase transition (during which protons and neutrons 
wcire assembled out of qiiarks and giuons), the dark matter freeze-o\it epoch 
(during which the dark matter decoupled from the cosmic plasma), neutrino 
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decoupling, electron-positron annihilation, and light-element nucleosynthesis 

(during which helium, deuterium and lithium were synthesized). The signa- 
tures that these processes left in the Universe can be used to constrain its 
parameters and answer the above questions. 

Half a million years after the big bang, hydrogen recombined and the 
Universe became transparent. The ultimate goal of observational cosmology 
is to image the entire history of the Universe since then. Currently, we have a 
snapshot of the Universe at recombination from the CMB, and detailed images 
of its evolution starting from an age of a billion years until the present time. 
The evolution between a million and a billion years has not been imaged as 
of yet. 

Within the next decade, NASA plans to launch an infrared space telescope 
(JWST) that will image the very first sources of light (stars and black holes) 
in the Universe, which are predicted theoretically to have formed in the first 
hundreds of millions of years. In parallel, there are several initiatives to con- 
struct large-aperture infrared telescopes on the ground with the same goal 
in mind"'^' The neutral hydrogen, relic from cosmological recombination, 
can be mapped in three-dimensions through its 21cm line even before the first 
galaxies formed [226] . Several groups are currently constructing low- frequency 
radio arrays in an attempt to map the initial inhomogcneitics as well as the 
process by which the hydrogen was re-ionized by the first galaxies. 

The next generation of ground-based telescopes will have a diameter of 
twenty to thirty meter. Together with JWST (that will not be affected by 
the atmospheric backgound) they will be able to image the first galaxies. 
Given that these galaxies also created the ionized bubbles around them, the 
same galaxy locations should correlate with bubbles in the neutral hydrogen 
(created by their UV emission). Within a decade it would be possible to 
explore the environmental influence of individual galaxies by using the two 
sets of instruments in concert [390] . 

The dark ingredients of the Universe can only be probed indirectly through 
a variety of luminous tracers. The distribution and nature of the dark matter 
are constrained by detailed X-ray and optical observations of galaxies and 
galaxy clusters. The evolution of the dark energy with cosmic time will be 
constrained over the coming decade by surveys of Type la supernovae, as well 
as surveys of X-ray clusters, up to a redshift of two. 

On large scales (> lOMpc) the power-spectrum of primordial density 
perturbations is already known from the measured microwave background 
anisotropics, galaxy surveys, weak Icnsing, and the Lya forest. Future pro- 
grams will refine current knowledge, and will search for additional trademarks 
of inflation, such as gravitational waves (through CMB polarization), small- 



^ http:/ /www. eso.org/projects/owl/ 
^ http://celt.ucolick.org/ 
^ http://www.gmto.org/ 
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Fig. 4. A sketch of the current design for the James Webb Space Tele- 
scope, the successor to the Hubble Space Telescope to be launched in 2011 
(http://www.jwst.nasa.gov/). The current design includes a primary mirror made 
of beryllium which is 6.5 meter in diameter as well as an instrument sensitivity that 
spans the full range of infrared wavelengths of 0.6-28/im and will allow detection 
of the first galaxies in the infant Universe. The telescope will orbit 1.5 million km 
from Earth at the Lagrange L2 point. 




Fig. 5. Artist conception of the design for one of the future giant telescopes that 
could probe the first generation of galaxies from the ground. The Giant Magellan 
Telescope (GMT) will contain seven mirrors (each 8.4 meter in diameter) and will 
have the resolving power equivalent to a 24.5 meter (80 foot) primary mirror. For 
more details see http://www.gmto.org/ 
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scale structure (through high-redshift galaxy surveys and 21cm studies), or 
the Gaussian statistics of the initial perturbations. 

The big bang is the only known event where particles with energies ap- 
proaching the Planck scale [(S-c^/G)^/^ ~ 10^^ GeV] interacted. It therefore 
offers prospects for probing the unification physics between quantum me- 
chanics and general relativity (to which string theory is the most-popular 
candidate). Unfortunately, the exponential expansion of the Universe during 
inflation erases memory of earlier cosmic epochs, such as the Planck time. 



3.3 Linear Gravitational Growth 



Observations of the CMB (e.g., Bennett et al. 1996 [36]) show that the Uni- 
verse at recombination was extremely uniform, but with spatial fluctuations 
in the energy density and gravitational potential of roughly one part in 10^. 
Such small fluctuations, generated in the early Universe, grow over time due 
to gravitational instability, and eventually lead to the formation of galaxies 
and the large-scale structure observed in the present Universe. 

As before, we distinguish between fixed and comoving coordinates. Using 
vector notation, the fixed coordinate r corresponds to a comoving position x = 
r/o. In a homogeneous Universe with density p, we describe the cosmological 
expansion in terms of an ideal pressureless fluid of particles each of which is 
at fixed x, expanding with the Hubble flow v = H{t)r where v = dr/dt. Onto 
this uniform expansion we impose small perturbations, given by a relative 
density perturbation 

^(x) = ^-l, (9) 

where the mean fluid density is p, with a corresponding peculiar velocity 
u = V — Hr. Then the fluid is described by the continuity and Euler equations 
in comoving coordinates [283, 284]: 

|^ + iv[(l+5)u]=0 (10) 

^ + ifu+-(u- V)u = — V(/). (11) 
at a a 

The potential cj) is given by the Poisson equation, in terms of the density 
perturbation: 

V^cj} = AnGpa^S . (12) 

This fluid description is valid for describing the evolution of collisionless 
cold dark matter particles until different particle streams cross. This "shell- 
crossing" typically occurs only after perturbations have grown to become non- 
linear, and at that point the individual particle trajectories must in general 
be followed. Similarly, baryons can be described as a pressureless fluid as long 
as their temperature is negligibly small, but non-linear collapse leads to the 
formation of shocks in the gas. 
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For small perturbations 6 <^ 1, the fluid equations can be linearized and 
combined to yield 

§+2Hf^=,.GpS. (13) 

This hncar equation has in general two independent solutions, only one of 
which grows with time. Starting with random initial conditions, this "growing 
mode" comes to dominate the density evolution. Thus, until it becomes non- 
linear, the density perturbation maintains its shape in comoving coordinates 
and grows in proportion to a growth factor D{t). The growth factor in the 
matter-dominated era is given by [283] 

^^^^ {^Aa^ + n^a + Hm)^^^ r a'^'^da' 

where we neglect fl^ when considering halos forming in the matter-dominated 
regime at z <C 10^. In the Einstein-de Sitter model (or, at high redshift, in 
other models as well) the growth factor is simply proportional to a(t). 

The spatial form of the initial density fluctuations can be described in 
Fourier space, in terms of Fourier components 

= y d^a;^(a;)e-*-^ . (15) 

Here we use the comoving wavcvcctor k, whose magnitude fc is the comov- 
ing wavenumber which is equal to 27t divided by the wavelength. The Fourier 
description is particularly simple for fluctuations generated by inflation (e.g., 
Kolb & Turner 1990 [205]). Inflation generates perturbations given by a Gaus- 
sian random field, in which different k-modes are statistically independent, 
each with a random phase. The statistical properties of the fiuctuations are 
determined by the variance of the different k-modes, and the variance is de- 
scribed in terms of the power spectrum P(fc) as follows: 

(4'5i:,) = (27r)^P(fc)j(3)(k-k') , (16) 

where J*-^-* is the three-dimensional Dirac delta function. The gravitational po- 
tential fluctuations are sourced by the density fluctuations through Poisson's 

equation. 

In standard models, inflation produces a primordial power-law spectrum 
P(fc) oc fc" with n ~ 1. Perturbation growth in the radiation-dominated and 
then matter-dominated Universe results in a modified final power spectrum, 
characterized by a turnover at a scale of order the horizon cH^^ at matter- 
radiation equality, and a small-scale asymptotic shape of P(fc) oc The 
overall amplitude of the power spectrum is not specified by current models of 
inflation, and it is usually set by comparing to the observed CMB temperature 
fluctuations or to local measures of large-scale structure. 
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Since density fluctuations may exist on all scales, in order to determine the 
formation of objects of a given size or mass it is useful to consider the statistical 
distribution of the smoothed density field. Using a window function W{r) 
normalized so that / d^rW{r) = 1, the smoothed density perturbation field, 
J d'^r5{x.)W{r), itself follows a Gaussian distribution with zero mean. For the 
particular choice of a spherical top-hat, in which W = Una, sphere of radius R 
and is zero outside, the smoothed perturbation field measures the fluctuations 
in the mass in spheres of radius R. The normalization of the present power 
spectrum is often specified by the value of as = cr{R = 8h~^Mpc). For the 
top-hat, the smoothed perturbation field is denoted Sn or Sm, where the mass 
M is related to the comoving radius _R by Af = AirpmR'^ /3, in terms of the 
current mean density of matter pm. The variance {Sm)^ is 

a\M)=a\R) = l^^^eP{k) 

where ji{x) = (sinx — .x cos 3;)/.t^ . The fimction (t{M) plays a crucial role in 
estimates of the abundance of collapsed objects, as we describe later. 

Species that decouple from the cosmic plasma (like the dark matter or the 
baryons) would show fossil evidence for acoustic oscillations in their power 
spectrum of inhomogeneities due to sound waves in the radiation fluid to 
which they were coupled at early times. This phenomenon can be understood 
as follows. Imagine a localized point-like perturbation from inflation at i = 0. 
The small perturbation in density or pressure will send out a sound wave that 
will reach the sound horizon c^t at any later time t. The perturbation will 
therefore correlate with its surroundings up to the sound horizon and all k- 
modes with wavelengths equal to this scale or its harmonics will be correlated. 
The scales of the perturbations that grow to become the first collapsed objects 
at z < 100 cross the horizon in the radiation dominated era after the dark 
matter decouples from the cosmic plasma. Next we consider the imprint of 
this decoupling on the smallest-scale structure of the dark matter. 

3.4 The Smallest-Scale Power Spectrum of Cold Dark Matter 

A broad range of observational data involving the dynamics of galaxies, the 
growth of large-scale structure, and the dynamics and nucleosynthesis of the 
Universe as a whole, indicate the existence of dark matter with a mean cosmic 
mass density that is ^ 5 times larger than the density of the baryonic matter 
[189, 348]. The data is consistent with a dark matter composed of weakly- 
interacting, massive particles, that decoupled early and adiabatically cooled 
to an extremely low temperature by the present time [189]. The Cold Dark 
Matter (CDM) has not been observed directly as of yet, although laboratory 
searches for particles from the dark halo of our own Milky- Way galaxy have 
been able to rc^strict the allowed parameter space for these particles. Since 
an alternative more-radical interpretation of the dark matter phenomenology 
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involves a modification of gravity [253] , it is of prime importance to find direct 
fingerprints of the CDM particles. One such fingerprint involves the small-scale 
structure in the Universe [158], on which we focus in this section. 

The most popular candidate for the CDM particle is a Weakly Interacting 
Massive Particle (WIMP). The lightest supcrsymmctric particle (LSP) could 
be a WIMP (for a review see [189]). The CDM particle mass depends on 
free parameters in the particle physics model but typical values cover a range 
around M ^100 GcV (up to values close to a TeV). In many cases the LSP 
hypothesis will be tested at the Large Hadron Collider (e.g. [33]) or in direct 
detection experiments (e.g. [16]). 

The properties of the CDM particles affect their response to the small- 
scale primordial inhomogeneities produced during cosmic inflation. The par- 
ticle cross-section for scattering off standard model fermions sets the epoch 
of their thermal and kinematic decoupling from the cosmic plasma (which 
is signiflcantly later than the time when their abundance freezes-out at a 
temperature T ~ M). Thermal decoupling is defined as the time when the 
temperature of the CDM stops following that of the cosmic plasma while 
kinematic decoupling is defined as the time when the bulk motion of the two 
species start to differ. For CDM the epochs of thermal and kinetic decou- 
pling coincide. They occur when the time it takes for collisions to change the 
momentum of the CDM particles equals the Hubble time. The particle mass 
determines the thermal spread in the speeds of CDM particles, which tends to 
smooth-out fluctuations on very small scales due to the free-streaming of parti- 
cles after kinematic decoupling [158, 159]. Viscosity has a similar effect before 
the CDM fluid decouples from the cosmic radiation fluid [182]. An important 
effect involves the memory the CDM fluid has of the acoustic oscillations of 
the cosmic radiation fluid out of which it decoupled. Here we consider the 
imprint of these acoustic oscillations on the small-scale power spectrum of 
density fluctuations in the Universe. Analogous imprints of acoustic oscilla- 
tions of the baryons were identified recently in maps of the CMB [348] , and the 
distribution of nearby galaxies [119]; these signatures appear on much larger 
scales, since the baryons decouple much later when the scale of the horizon is 
larger. The discussion in this section follows Loeb & Zaldarriaga (2005) [228]. 
Formalism 

Kinematic decoupling of CDM occurs during the radiation-dominated era. 
For example, if the CDM is made of neutralinos with a particle mass of ~ 
100 GeV, then kinematic decoupling occurs at a cosmic temperature of Ta ~ 
10 MeV [182, 87]. As long as T<j < 100 MeV, we may ignore the imprint of 
the QCD phase transition (which transformed the cosmic quark-gluon soup 
into protons and neutrons) on the CDM power spectrum [321]. Over a short 
period of time during this transition, the pressure does not depend on density 
and the sound speed of the plasma vanishes, resulting in a significant growth 
for perturbations with periods shorter than the length of time over which 
the sound speed vanishes. The transition occurs when the temperature of 
the cosmic plasma is ~ 100 — 200 MeV and lasts for a small fraction of the 
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Hubble time. As a result, the induced modifications are on scales smaller than 
those we arc considering here and the imprint of the QCD phase transition is 
washed-out by the effects we calculate. 

At early times the contribution of the dark matter to the energy density 
is negligible. Only at relatively late times when the cosmic temperature drops 
to values as low as ~ 1 eV, matter and radiation have comparable energy 
densities. As a result, the dynamics of the plasma at earlier times is virtually 
unaffected by the presence of the dark matter particles. In this limit, the 
dynamics of the radiation determines the gravitational potential and the dark 
matter just responds to that potential. We will use this simplification to obtain 
analytic estimates for the behavior of the dark matter transfer function. 

The primordial inflationary fluctuations lead to acoustic modes in the radi- 
ation fluid during this era. The interaction rate of the particles in the plasma is 
so high that we can consider the plasma as a perfect fluid down to a comoving 
scale, 

A/~77d/\/A^ ; N^natd, (18) 

where rjd = Jq'' dt/a{t) is the conformal time (i.e. the comoving size of the 
horizon) at the time of CDM decoupling, t^; cr is the scattering cross section 
and n is the relevant particle density. (Throughout this section we set the 
speed of light and Planck's constant to unity for brevity.) The damping scale 
depends on the species being considered and its contribution to the energy 
density, and is the largest for neutrinos which are only coupled through weak 
interactions. In that case N ~ {T/T^Y where ~ 1 MeV is the temperature 
of neutrino decoupling. At the time of CDM decoupling N ^ M/Td ~ 10^ for 
the rest of the plasma, where M is the mass of the CDM particle. Here we 
will consider modes of wavelength larger than A / , and so we neglect the effect 
of radiation diffusion damping and treat the plasma (without the CDM) as a 
perfect fluid. 

The equations of motion for a perfect fluid during the radiation era can 
be solved analytically. We will use that solution here, following the notation 
of Dodelson [109]. As usual we Fourier decompose fluctuations and study 
the behavior of each Fourier component separately. For a mode of comoving 
wavenumber k in Newtonian gauge, the gravitational potential fluctuations 
are given by: 

- '"^4 (a;^)3 

where u ^ k/ \/3 is the frequency of a mode and #p is its primordial ampli- 
tude in the limit ?y — » 0. In this section we use conformal time r] = J dt/a{t) 
with a{t) DC t^/^ during the radiation-dominated era. Expanding the temper- 
ature anisotropy in multipole moments and using the Boltzmann equation 
to describe their evolution, the monopole 0o and dipole 0i of the photon 
distribution can be written in terms of the gravitational potential as [109]: 



(19) 
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e, = 4(4>. + i*) (20) 

where x = krj and a prime denotes a derivative with respect to x. 

The solutions in equations (19) and (20) assume that both the sound speed 
and the number of relativistic degrees of freedom arc constant over time. As 
a result of the QCD phase transition and of various particles becoming non- 
relativistic, both of these assumptions are not strictly correct. The vanishing 
sound speed during the QCD phase transition provides the most dramatic 
effect, but its imprint is on scales smaller than the ones we consider here 
because the transition occurs at a significantly higher temperature and only 
lasts for a fraction of the Hubble time [321]. 

Before the dark matter decouples kinematically, we will treat it as a fluid 
which can exchange momentum with the plasma through particle collisions. 
At early times, the CDM fluid follows the motion of the plasma and is involved 
in its acoustic oscillations. The continuity and momentum equations for the 
CDM can be written as: 

Sc + 9c = 

Oc + -9c = k'^c% - k'^cJc - + T-\ei - 9c) (21) 
a 

where a dot denotes an r;-derivative, 5c is the dark matter density perturba- 
tion, 9c is the divergence of the dark matter velocity field and ctc denotes the 
anisotropic stress. In writing these equations we have followed Rcf. [230]. The 
term t^^(6'i — Oc) encodes the transfer of momentun between the radiation 
and CDM fluids and provides the collisional rate of momentum transfer, 

= na^a, (22) 

with n being the number density of particles with which the dark matter is 
interacting, cr(T) the average cross section for interaction and M the mass of 
the dark matter particle. The relevant scattering partners are the standard 
model leptons which have thermal abundances. For detailed expressions of 
the cross section in the case of supersymmetric (SUSY) dark matter, see Refs. 
[87, 159]. For our purpose, it is sufficient to specify the typical size of the cross 
section and its scaling with cosmic time, 

where the coupling mass M„ is of the order of the weak-interaction scale 
(~ 100 GeV) for SUSY dark matter. This equation should be taken as the 
definition of , as it encodes all the imccrtaintics in the details of the par- 
ticle physics model into a single parameter. The temperature dependance of 
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the averaged cross section is a result of the available phase space. Our re- 
sults arc quite insensitive to the details other than through the decoupling 
time. Equating r~^/a to the Hubble expansion rate gives the temperature of 
kinematic decoupling: 

f M'^mV^ ( M„ \ f M V'"^ 

''"{-Jur) =='°"'nio«(iev)(iiM5;v) ■ f^"' 

The term k^cj.6c in Eq. (21) results from the pressure gradient force and Cs 
is the dark matter sound speed. In the tight coupling limit, Tc <C H^^ we find 
that w fcT/M and that the shear term is fc^CTc « fvC^TcOc- Here /„ and fc 
are constant factors of order unity. We will find that both these terms make a 
small difference on the scales of interest, so their precise value is unimportant. 

By combining both equations in (21) into a single equation for Sc we get 

= S{x)-3F^{x)<P' + ^{W',-6'J, (25) 

where Xd = kr]d and rjd denotes the time of kinematic decoupling which can 
be expressed in terms of the decoupling temperature as, 

oc M-^M-i/'', (26) 

with To = 2.7K being the present-day CMB temperature and Zd being the 
redshift at kinematic decoupling. We have also introduced the source function, 

S(x) = -3^" + # - (27) 

X 

For X <^ Xd, the dark matter sound speed is given by 

clix)=clixd)^, (28) 

X 

where cl{xd) is the dark matter sound speed at kinematic decoupling (in units 
of the speed of light). 

In writing (28) we have assumed that prior to decoupling the temperature of 
the dark matter follows that of the plasma. For the viscosity term we have, 

F,{x) = Ul{xd)xl(^^y . (30) 
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Free streaming after kinematic decoupling 

In the limit of the collision rate being much slower than the Hubble expan- 
sion, the CDM is decoupled and the evolution of its perturbations is obtained 
by solving a Boltzman equation: 



d£ _^ dndf_ ^ dqj df 

dr] dr] dvi drj dqi ' 



(31) 



where f{r,q,r]) is the distribution function which depends on position, co- 
moving momentum q, and time. The comoving momentum 3-components are 
dxi /dr] = qi /a. We use the Boltzman equation to find the evolution of modes 
that are well inside the horizon with a; » 1. In the radiation era, the grav- 
itational potential decays after horizon crossing (sec Eq. 19). In this limit 
the comoving momentum remains constant, dqi/drj = and the Boltzman 
equation becomes, 

drj a dvi 

We consider a single Fourier mode and write / as, 



= 0. 



(32) 



where fo{Q) is the unperturbed distribution. 



fo{Q) 



( M \ 



3/2 



riCDM exp 

V ^"^J CDM / 



1 Mq^ 



2 TcDM 



(33) 



(34) 



where ncoM and Tcdm are the present-day density and temperature of the 
dark matter. 

Our approach is to solve the Boltzman equation with initial conditions 
given by the fluid solution at a time 77* (which will depend on k) . The simplified 
Boltzman equation can be easily solved to give SF{q,ri) as a function of the 
initial conditions (5i?(q, 77,), 



Spiq, v) = ^F{q, V*) exp[-iqr • fe 



V* 



(35) 



The CDM overdensity Sc can then be expressed in terms of the perturba- 
tion in the distribution function as, 



1 



d^q foiq) ^F{q,v)- 



ncDM 

We can use (35) to obtain the evolution of 5c in terms of its value at 77,, 

IP 
V 



(36) 



5c{ri) = exp 



1 2/ »7 ^ 



CI I T f V ^ 

'5k + ^k^*ln(-) 



(37) 
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where kj"^ = (Td/M)rid- The exponential term is responsible for the damp- 
ing of perturbations as a result of free streaming and the disperwion of the 
CDM particles after they decouple from the plasma. The above expression is 
only valid during the radiation era. The free streaming scale is simply given 
by j dt{v / a) cx J dta~^ which grows logarithmically during the radiation era 
as in equation (37) but stops growing in the matter era when a oc t"^^^. 

Equation (37) can be used to show that even during the free streaming 
epoch, Sc satisfies equation (25) but with a modified sound speed and viscous 
term. For x ^ Xd one should use, 



l + xlclixd)!!^^- . 

Xd 



F,ix) = 2cl{xd)xlln(^^) (38) 



The differences between the above scalings and those during the tight coupling 
regime are a result of the fact that the dark matter temperature stops follow- 
ing the plasma temperature but rather scales as after thermal decoupling, 
which coincides with the kinematic decoupling. We ignore the efi'ects of heat 
transfer during the fluid stage of the CDM because its temperature is con- 
trolled by the much larger heat reservoir of the radiation-dominated plasma 
at that stage. 

To obtain the transfer function we solve the dark matter fluid equation 
until decoupling and then evolve the overdensity using equation (37) up to the 
time of matter-radiation equality. In practice, we use the fluid equations up to 
a;* = 10 max{xd, 10) so as to switch into the free streaming solution well after 
the gravitational potential has decayed. In the fluid equations, we smoothly 
match the sound speed and viscosity terms at x = a;^. As mentioned earlier, 
because Cs{xd) is so small and we are interested in modes that are comparable 
to the size of the horizon at decoupling, i.e. Xd ~ few, both the dark matter 
sound speed and the associated viscosity play only a minor role, and our 
simplifled treatment is adequate. 

In Figure 6 we illustrate the time evolution of modes during decoupling 
for a variety of k values. The situation is clear. Modes that enter the horizon 
before kinematic decoupling oscillate with the radiation fluid. This behavior 
has two important effects. In the absence of the coupling, modes receive a 
"kick" by the source term S{x) as they cross the horizon. After that they 
grow logarithmically. In our case, modes that entered the horizon before kine- 
matic decoupling follow the plasma oscillations and thus miss out on both the 
horizon "kick" and the beginning of the logarithmic growth. Second, the de- 
coupling from the radiation fluid is not instantaneous and this acts to further 
damp the amplitude of modes with Xd ^ 1. This effect can be understood as 
follows. Once the oscillation frequency of the mode becomes high compared 
to the scattering rate, the coupling to the plasma effectively damps the mode. 
In that limit one can replace the forcing term 6*0 by its average value, which 
is close to zero. Thus in this regime, the scattering is forcing the amplitude 
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Fig. 6. The normahzed ampHtude of CDM fluctuations 5/<?>p for a variety of modes 
with comoving wavenumbers log{krjd) — (0,1/3,2/3, 1,4/3,5/3,2) as a function of 
X = krj, where rj = j^dt/ait) is the conformal time coordinate. The dashed line 
shows the temperature monopole Wo and the uppermost (dotted) curve shows the 
evolution of a mode that is uncoupled to the cosmic plasma. 

of the dark matter oscillations to zero. After kinematic decoupling the modes 
again grow logarithmicaUy but from a very reduced amplitude. The coupling 
with the plasma induces both oscillations and damping of modes that entered 
the horizon before kinematic decoupling. This damping is different from the 
free streaming damping that occurs after kinematic decoupling. 

In Figure 7 we show the resulting transfer function of the CDM overden- 
sity. The transfer function is defined as the ratio between the CDM density 
perturbation amplitude 5c when the effect of the coupling to the plasma is 
included and the same quantity in a model where the CDM is a perfect fluid 
down to arbitrarily small scales (thus, the power spectrum is obtained by 
multiplying the standard result by the square of the transfer function). This 
function shows both the oscillations and the damping signature mentioned 
above. The peaks occur at multipoles of the horizon scale at decoupling, 

kpeak = (8, 15.7, 24.7, ..)7^2^ (X (39) 

This same scale determines the "oscillation" damping. The free streaming 
damping scale is, 

VdCdiVd) ln{rieq/r]d) cx ^' \n{Td/Tcq), (40) 

ToT^' 

where T^q is the temperature at matter radiation equality, Toq « 1 eV. The 
free streaming scale is parametrically different from the "oscillation" damping 
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Fig. 7. Transfer function of the CDM density perturbation amplitude (normalized 
by the primordial amplitude from inflation). We show two cases: (i) Td/M — 10~* 
and Td/Toq = 10^; (ii) Td/M = 10'^ and Td/Teq = 10^. In each case the oscillatory 
curve is our result and the other curve is the free-streaming only result that was 
derived previously in the literature [4,7,8]. 



scale. However for our fiducial choice of parameters for the CDM particle they 
roughly coincide. 

The CDM damping scale is significantly smaller than the scales observed 
directly in the Cosmic Microwave Background or through large scale structure 
surveys. For example, the ratio of the damping scale to the scale that entered 
the horizon at Matter-radiation equality is rjd/rjeq ~ Teq/Td ~ 10~^ and 
to our present horizon rjd/rjo ^ (T'eqTo)^/^/Td ~ 10^^. In the context of 
infiation, these scales were created 16 and 20 e-folds apart. Given the large 
extrapolation, one could certainly imagine that a change in the spectrum could 
alter the shape of the power spectrum around the damping scale. However, for 
smooth inflaton potentials with small departures from scale invariance this is 
not likely to be the case. On scales much smaller than the horizon at matter 
radiation equality, the spectrum of perturbations density before the effects of 
the damping are included is approximately. 



{n - 1) ln(A;77e,) + ln(fc77e,)^ + 



A (k) cx exp 

xln"(fc77e,/8) (41) 



where the first term encodes the shape of the primordial spectrum and the 
second the transfer function. Primordial departures from scale invariance are 
encoded in the slope n and its running a. The effective slope at scale k is 
then, 

^ = (n-l)+oln(fcW + ^^- (42) 
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For typical values of (n — 1) ~ 1/60 and a ~ 1/60^ the slope is still positive 
at A: ^ 77^^^, so the cut-off in the power will come from the effects we calcu- 
late rather than from the shape of the primordial spectrum. However given 
the large extrapolation in scale, one should keep in mind the possibility of 
significant effects resulting from the mechanisms that generates the density 
perturbations. 

Implications We have found that acoustic oscillations, a relic from the epoch 

when the dark matter coupled to the cosmic radiation fluid, truncate the CDM 
power spectrum on a comoving scale larger than effects considered before, such 
as free-streaming and viscosity [158, 159, 182]. For SUSY dark matter, the 
minimum mass of dark matter clumps that form in the Universe is therefore 
increased by more than an order of magnitude to a value of 



where Pcrit = {Hq/SttG) = 9 x 10""^" g cm~'^ is the critical density today, and 
i?M is the matter density for the concordance cosmological model [348]. We 

define the cut-off wavenumber fccut as the point where the transfer function 
first drops to a fraction 1/e of its value at k ^ 0. This corresponds to kcut ~ 



Recent numerical simulations [105, 146] of the earliest and smallest objects 
to have formed in the Universe (see Fig. 3.4), need to be redone for the modi- 
fied power spectrum that we calculated in this section. Although it is difficult 
to forecast the effects of the acoustic oscillations through the standard Press- 
Schechter formalism [291], it is likely that the results of such simulations will 
be qualitatively the same as before except that the smallest clumps would 
have a mass larger than before (as given by Eq. 43). 

Potentially, there are several observational signatures of the smallest CDM 
clumps. As pointed out in the literature [105, 353], the smallest CDM clumps 
coiild produce 7-rays through dark-matter annihilation in their inner density 
cusps, with a flux in excess of that from nearby dwarf galaxies. If a sub- 
stantial fraction of the Milky Way halo is composed of CDM clumps with 
a mass ~ 10~^Mq, the nearest clump is expected to be at a distance of 
^ 4 X 10^^ cm. Given that the characteristic speed of such clumps is a few 
hundred km s^^, the 7-ray flux would therefore show temporal variations on 
the relatively long timescale of a thousand years. Passage of clumps through 
the solar system should also induce fluctuations in the detection rate of CDM 
particles in direct search experiments. Other observational effects have rather 



* Our definition of the cut-off mass follows the convention of the Jeans mass, which 
is defined as the mass enclosed within a sphere of radius Aj/2 where Aj = 27r/fej 
is the Jeans wavelength [168]. 




(43) 



3.3 rj^\ 
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Fig. 8. A slice through a numerical simulation of the first dark matter condensations 
to form in the Universe. Colors represent the dark matter density at z — 26. The 
simulated volume is 60 comoving pc on a side, simulated with 64 million particles 
each weighing 1.2 x 10"^°Mq (!). (from Diemand, Moore, & Stadel 2005 [105]). 

limited prospects for detectability. Because of their relatively low-mass and 
large size (~ 10" cm), the CDM clump s are too diffuse to produce any gravita- 
tional lensing signatures (including /emto-lensing [161]), even at cosmological 
distances. 

The smallest CDM clumps should not affect the intergalactic baryons 
which have a much larger Jeans mass. However, once objects above ~ IO^Mq 
start to collapse at redshifts z < 30, the baryons would be able to cool inside 
of them via molecular hydrogen transitions and the interior baryonic Jeans 
mass would drop. The existence of dark matter clumps could then seed the 
formation of the first stars inside these objects [66]. 

3.5 Structure of the Baryons 

Early Evolution of Baryonic Perturbations on Large Scales 

The baryons are coupled through Thomson scattering to the radiation 
fluid until they become neutral and decouple. After cosmic recombination, 
they start to fall into the potential wells of the dark matter and their early 
evolution was derived by Barkana & Loeb (2005) [29]. 

On large scales, the dark matter (dm) and the baryons (b) are affected only 
by their combined gravity and gas pressure can be ignored. The evolution of 
sub-horizon linear perturbations is described in the matter-dominated regime 
by two coupled second-order differential equations [284] : 

6dm + 2ff(5dm = 4:TTGp,n (/b^b + /dm'^dm) , 

+ 2HSb = 4TTGp,n ifbSb + fdnM , (44) 
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where ^dm(i) and Sh{t) are the perturbations in the dark matter and baryons, 
respectively, the derivatives arc witli respect to cosmic time t, H{t) = a/a is 
the Hubble constant with a = {1 + z)~^ , and we assume that the mean mass 
density pm{t) is made up of respective mass fractions /dm and /b = 1 — /dm- 
Since these linear equations contain no spatial gradients, they can be solved 
spatially for 5dm (x, t) and (5b (x, t) or in Fourier space for 5dm (k, t) and 5b (k, t). 
Defining 5tot = fh^h + /dm^dm and 5b- = 5b - 5tot , we find 

5tot + 2ff 5tot = A-nGpm^tot , 

5b- + 2if5b- = . (45) 

Each of these equations has two independent solutions. The equation for 5tot 
has the usual growing and decaying solutions, which we denote Di{t) and 
D4{t), respectively, while the 5b- equation has solutions D2{t) and D:i{t)-, we 
number the solutions in order of declining growth rate (or increasing decay 
rate). We assume an Einstein-de Sitter, matter- dominated Universe in the red- 
shift range z = 20 150, since the radiation contributes less than a few percent 
at z < 150, while the cosmological constant and the curvature contribute to 
the energy density less than a few percent at ^ > 3. In this regime a oc t^^^ and 
the solutions are Di(t) = a/ui and D4^{t) = (a/ai)^^^^ for 5tot, and D2{t) = 1 
and D3{t) = {a/ai)~^^'^ for 5b-, where we have normalized each solution to 
unity at the starting scale factor Oj, which we set at a redshift Zi = 150. The 
observable baryon perturbation can then be written as 

4 

5b (k, t) = 5b- + 5tot = 13 (k) Dm {t) , (46) 

m=l 

and similarly for the dark matter perturbation, 

4 

5dm = (Stot - fbSb) = J2 ^^C*) ' (4^) 

where Cj = Di for i = 1,4 and Ci = — (/b//dm)-Di for i = 2,3. We may 
establish the values of 5„(k) by inverting the 4x4 matrix A that relates the 
4- vector (5i, 52, 53, 54) to the 4- vector that represents the initial conditions 

(5b, 5dm, 5b, 5dm) at the initial time. 

Next we describe the fluctuations in the sound speed of the cosmic gas 
caused by Compton heating of the gas, which is due to scattering of the 
residual electrons with the CMB photons. The evolution of the temperature 
T of a gas element of density pj, is given by the first law of thermodynamics: 

3 

dQ = - kdT - kTd log pb , (48) 
where dQ is the heating rate per particle. Before the first galaxies formed, 



First Light 23 



10^ p 1 r I r — I — I 1 1 1 1 I I I — I — I — I 1 1 




l + z 

Fig. 9. Redshift evolution of the amphtudes of the independent modes of the den- 
sity perturbations (upper panel) and the temperature perturbations (lower panel), 
starting at redshift 150 (from Barkana & Loeb 2005 [29]). We show m — 1 (solid 
curves), m = 2 (short-dashed curves), m = 3 (long-dashed curves), m = 4 (dotted 
curves), and m = (dot-dashed curve). Note that the lower panel shows one plus 
the mode amplitude. 



^ = 4^fc(T^_T)p^^,(i) , (49) 

at TUe 

where ax is the Thomson cross-section, Xe(t) is the electron fraction out of 
the total number density of gas particles, and is the CMB energy density 
at a temperature T^. In the redshift range of interest, we assume that the 
photon temperature (Tj = T^/a) is spatially uniform, since the high sound 
speed of the photons (i.e., c/\/3) suppresses fluctuations on the sub-horizon 
scales that we consider, and the horizon-scale ^ 10~^ fluctuations imprinted 
at cosmic recombination are also negligible compared to the smallWe estab- 
lish the values of ^^(k) by inverting the 4x4 matrix A that relates the 
4- vector (5i, ^2, ^3, ^4) to the 4- vector that represents the initial conditions 

((^b, <5dm, (^b, <5dm) at the initial time, er-scale fluctuations in the gas density 
and temperature. Fluctuations in the residual electron fraction x^it) are even 
smaller. Thus, 

dT 2 dlogpb ^ Xe{t) _4 
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Fig. 10. Power spectra and initial perturbation amplitudes versus wavenumber 
(from [29]). The upper panel shows Sh (solid curves) and 5dm (dashed curves) at 
z = 150 and 20 (from bottom to top). The lower panel shows the initial {z — 150) 
amplitudes of Si (solid curve), S2 (short-dashed curve), S3 (long-dashed curve), 84 
(dotted curve), and Sxiti) (dot-dashed curve). Note that if Si is positive then so are 
S3 and Sriti), while ^2 is negative at all k, and ^4 is negative at the lowest k but is 
positive at fc > 0.017 Mpc"^ 



where = p°(8crT c/3me) = 8.55 x 10 ^'^yr ^. After cosmic recombination, 
Xe{t) changes due to the slow recombination rate of the residual ions: 

^^^^aB{T)xl{t)nH{l+y) , (51) 

where asiT) is the case-B recombination coefficient of hydrogen, nn is the 
mean number density of hydrogen at time t, and y — 0.079 is the helium to 
hydrogen number density ratio. This yields the evolution of the mean tem- 
perature, df/dt = -2Hf + Xe{t)t-'^ {T^-f)a-'^. In prior analyses [284, 230] 
a spatially uniform speed of sound was assumed for the gas at each redshift. 
Note that we refer to 5p/dp as the square of the sound speed of the fluid, 
where 5p is the pressure perturbation, although we are analyzing perturba- 
tions driven by gravity rather than sound waves driven by pressure gradients. 

Instead of assuming a uniform sound speed, we find the first-order pertur- 
bation equation. 



ddr 2 dSb Xe{t)T^ _4 
dt 3 dt T 



St , (52) 
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where we defined the fractional temperature perturbation St- Like the density 
perturbation equations, this equation can be solved separately at each x or 
at each k. Furthermore, the solution Sxit) is a linear functional of Si,{t) [for a 
fixed function Xe(t)]. Thus, if we choose an initial time ti then using Eq. (46) 
we can write the solution in Fourier space as 

4 

~St{K = E ^™(*) + U) Dl{t) , (53) 

m— 1 

where D^{t) is the solution of Eq. (52) with (5^ = at ti and with the 
perturbation mode Dm{t) substituted for (5f,(i), while Dj^^t) is the solution 
with no perturbation Sbit) and with (5t = 1 at t^. By modifying the CMBFAST 
code (http://www.cmbfast.org/), we can numerically solve Eq. (52) along with 
the density perturbation equations for each k down to Zi — 150, and then 
match the solution to the form of Eq. (53). 
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Fig. 11. Relative sensitivity of perturbation amplitudes at z = 150 to cos- 
mological parameters (from [29]). For variations in a parameter x, we show 
dlog-^P(fc)/dlog(a;). We consider variations in Odmh^ (upper panel), in Obh^ (mid- 
dle panel), and in the Hubble constant h (lower panel). When we vary each param- 
eter we fix the other two, and the variations are all carried out in a flat i^totai — 1 
universe. We show the sensitivity of 5i (solid curves), 62 (short-dashed curves), S3 
(long-dashed curves), and Sxiti) (dot-dashed curves). 
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Figure 9 shows the time evolution of the various independent modes that 

make up the perturbations of density and temperature, starting at the time 
ti corresponding to Zj = 150. is identically zero since D2{t) = 1 is con- 

stant, while D'^{t) and Dj{t) are negative. Figure 10 shows the amplitudes 
of the various components of the initial perturbations. We consider comov- 
ing wavevectors k in the range 0.01 - 40 Mpc~^, where the lower limit is set 
by considering sub-horizon scales at z = 150 for which photon perturbations 
are negligible compared to 5dm and Sh, and the upper limit is set by requir- 
ing baryonic pressure to be negligible compared to gravity. 62 and ^3 clearly 
show a strong signature of the large-scale baryonic oscillations, left over from 
the era of the photon-baryon fluid before recombination, while Si, S4, and St 
carry only a weak sign of the oscillations. For each quantity, the plot shows 
[fc^P(fc)/(27r^)]^/^, where P{k) is the corresponding power spectrum of fluc- 
tuations. Si is already a very small correction at z = 150 and declines quickly 
at lower redshift, but the other three modes all contribute significantly to S^, 
and the Sriti) term remains significant in Srit) even at z < 100. Note that 
at z = 150 the temperature perturbation St has a different shape with re- 
spect to k than the baryon perturbation S^,, showing that their ratio cannot 
be described by a scale-independent speed of sound. 

The power spectra of the various perturbation modes and of Sxiti) depend 
on the initial power spectrum of density fluctuations from inflation and on the 
values of the fundamental cosmological parameters (I7dm, ^b, ^a, and h). If 
these independent power spectra can be measured through 21cm fluctuations, 
this will probe the basic cosmological parameters through multiple combina- 
tions, allowing consistency checks that can be used to verify the adiabatic 
nature and the expected history of the perturbations. Figure 11 illustrates 
the relative sensitivity of a/ P{k) to variations in f^dmh'^, fii,h^, and h, for the 
quantities 5i, (52, (^3, and 5T{ti)- Not shown is ^4, which although it is more 
sensitive (changing by order imity due to 10% variations in the parameters), 
its magnitude always remains much smaller than the other modes, making it 
much harder to detect. Note that although the angular scale of the baryon 
oscillations constrains also the history of dark energy through the angular 
diameter distance, we have focused here on other cosmological parameters, 
since the contribution of dark energy relative to matter becomes negligible at 
high redshift. 

Cosmological Jeans Mass 

The Jeans length Aj was originally defined (Jeans 1928 [187]) in Newtonian 

gravity as the critical wavelength that separates oscillatory and exponentially- 
growing density perturbations in an infinite, uniform, and stationary distri- 
bution of gas. On scales (. smaller than Aj, the sound crossing time, (./cg is 
shorter than the gravitational free-fall time, {Gp)~^/'^, allowing the build-up 
of a pressure force that counteracts gravity. On larger scales, the pressure 
gradient force is too slow to react to a build-up of the attractive gravita- 
tional force. The Jeans mass is defined as the mass within a sphere of radius 
Aj/2, Mj = (47r/3)/9(Aj/2)^. In a perturbation with a mass greater than Mj, 
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the self-gravity cannot be supported by the pressure gradient, and so the gas 
is unstable to gravitational collapse. The Newtonian derivation of the Jeans 
instability suffers from a conceptual inconsistency, as the unperturbed gravi- 
tational force of the uniform background must induce bulk motions (compare 
to Binncy & Tremaine 1987 [43]). However, this inconsistency is remedied 
when the analysis is done in an expanding Universe. 

The perturbative derivation of the Jeans instability criterion can be car- 
ried out in a cosmological setting by considering a sinusoidal perturbation 
superposed on a uniformly expanding background. Here, as in the Newtonian 
limit, there is a critical wavelength Aj that separates oscillatory and growing 
modes. Although the expansion of the background slows down the exponential 
growth of the amplitude to a power-law growth, the fundamental concept of 
a minimum mass that can collapse at any given time remains the same (see, 
e.g. Kolb & Turner 1990 [205]; Peebles 1993 [284]). 

We consider a mixture of dark matter and baryons with density parameters 
— Pdm/Pc and = p\,/ pc, where pdm is the average dark matter density, 
Pb is the average baryonic density, pc is the critical density, and f^l^ + f2^ = 
is given by equation(83). We also assume spatial fluctuations in the gas 
and dark matter densities with the form of a single spherical Fourier mode on 
a scale much smaller than the horizon, 

Pdrajr, t)- Pdmjt) . sin(fcr) 

— J^) — -''-^'^^^ ^''^ 

p,{r,t) - p,{t) ^ sjn^ ^ 
ph(t) kr 

where Pdm(i) and Ph{t) are the background densities of the dark matter and 
baryons, (5dm(i) and 6h{t) are the dark matter and baryon overdensity ampli- 
tudes, r is the comoving radial coordinate, and k is the comoving perturbation 
wavenumber. We adopt an ideal gas equation-of-state for the baryons with a 
specific heat ratio 7=5/3. Initially, at time t = ti, the gas temperature is uni- 
form Tb(r, fi)=Ti, and the perturbation amplitudes are small (5dm,i,(5b,i -C 1. 
We define the region inside the first zero of sin(fcr)/(fcr), namely < kr < n, 
as the collapsing "object" . 

The evolution of the temperature of the baryons Tb (r, t) in the lin- 
ear regime is determined by the coupling of their free electrons to the 
CMB through Compton scattering, and by the adiabatic expansion of the 
gas. Hence, Tb(r, t) is generally somewhere between the CMB temperature, 
T-y oc (1 -I- z)~^ and the adiabatically-scaled temperature Tad oc (1 + z)~^. In 
the limit of tight coupling to T^, the gas temperature remains uniform. On 
the other hand, in the adiabatic limit, the temperature develops a gradient 
according to the relation 

Tb oc p'^-'\ (56) 

The evolution of a cold dark matter overdensity, Sdmit), in the linear regime 
is described by the equation (44), 
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5dm + 2if(5dm = {^h^h + ^^dm^dm) (57) 

whereas the evolution of the overdensity of the baryons, 5]^{t), with the inclu- 
sion of their pressure force is described by (see §9.3.2 of [205]), 

5b+2HS,, = -H^ni,5b + Odra6dm) -{-) (-) Ub + -/3[(5b-M 

2 iirup \a J \ a / \ 3 

(58) 

Here, H{t) = a/a is the Hubble parameter at a cosmological time t, and 
fj, = 1.22 is the mean molecular weight of the neutral primordial gas in atomic 
units. The parameter /3 distinguishes between the two limits for the evolution 
of the gas temperature. In the adiabatic limit (3 = 1, and when the baryon 
temperature is uniform and locked to the background radiation, /3 = 0. The 
last term on the right hand side (in square brackets) takes into account the 
extra pressure gradient force in \7{py,T) = (TVpb + PbVT), arising from the 
temperature gradient which develops in the adiabatic limit. The Jeans wave- 
length Aj = 27r/fcj is obtained by setting the right-hand side of equation (58) 
to zero, and solving for the critical wavenumber fcj. As can be seen from equa- 
tion (58), the critical wavelength Aj (and therefore the mass Mj) is in general 
time-dependent. We infer from equation (58) that as time proceeds, pertur- 
bations with increasingly smaller initial wavelengths stop oscillating and start 
to grow. 

To estimate the Jeans wavelength, we equate the right-hand-side of equa- 
tion (58) to zero. Wc further approximate Sb ~ <5dm, and consider suffi- 
ciently high rcdshifts at which the Universe is matter dominated and flat, 
(l+z) > max[(l-l2™-l?yi)/^2m, (^^A/^^m)^/^]. Inthis regime, i^b ~ 1^ 

H « 2/(3t), and a = (l + z)~^ ~ (3i^o\/^t^/2)^/^^^/^ where f2m = i7dm + i^6 
is the total matter density parameter. Following cosmological recombination 
at z w lO'', the residual ionization of the cosmic gas keeps its temperature 
locked to the CMB temperature (via Compton scattering) down to a redshift 
of [284] 

1 + w 160(j76/iV0.022)2/5 . (59) 
In the redshift range between recombination and zt, /? = and 

fcj = (27r/Aj) = [2kT^{0)/3ijmp]-^/^^/Q^Ho , (60) 

so that the Jeans mass is therefore redshift independent and obtains the value 
(for the total mass of baryons and dark matter) 

Based on the similarity of Mj to the mass of a globular cluster, Peebles & 
Dicke (1968) [281] suggested that globular clusters form as the first generation 
of baryonic objects shortly after cosmological recombination. Peebles & Dicke 
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assumed a baryonic Universe, with a nonlinear fluctuation aniphtude on sniaU 
scales at z ~ 10'^, a model which has by now been ruled out. The lack of a 
dominant mass of dark matter inside globular clusters makes it unlikely that 
they formed through direct cosmological collapse, and more likely that they 
resulted from fragmentation during the process of galaxy formation. 




1 10 100 1000 

(1 + z) 

Fig. 12. Thermal history of the baryons, left over from the big bang, before the first 
galaxies formed. The residual fraction of free electrons couple the gas temperture 
Tgas to the cosmic microwave background temperature [T^, cx (1 + ^)] until a redshift 
z ^ 200. Subsequently the gas temperature cools adiabatically at a faster rate 
[Tgas oc (1 + z)^\. Also shown is the spin temperature of the 21cm transition of 
hydrogen Ts which interpolates between the gas and radiation temperature and will 
be discussed in detail later in this review. 



At z < zt, the gas temperature declines adiabatically as [(1 + z)l(\ + Zf)]^ 
(i.e., /3 = 1) and the total Jeans mass obtains the value, 

Afj ^ 4.54 X 103 -ili- — - Mo. (62) 



^ 0.15 j VO-022 ) V 10 , 

It is not clear how the value of the Jeans mass derived above relates to the 
mass of collapsed, bound objects. The above analysis is perturbative (Eqs. 57 
and 58 are valid only as long as i5b and 5dm are much smaller than unity) , and 
thus can only describe the initial phase of the collapse. As (5b and (5dm grow and 
become larger than unity, the density profiles start to evolve and dark matter 
shells may cross baryonic shells [167] due to their different dynamics. Hence 
the amount of mass enclosed within a given baryonic shell may increase with 
time, until eventually the dark matter pulls the baryons with it and causes 
their collapse even for objects below the Jeans mass. 

Even within linear theory, the Jeans mass is related only to the evolution 
of perturbations at a given time. When the Jeans mass itself varies with time. 
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the overall suppression of the growth of perturbations depends on a time- 
weighted Jeans mass. Gncdin & Hui (1998) [150] showed that the correct 
time-weighted mass is the filtering mass Mp = (47r/3) p{2TTa/kF)^ , in terms 
of the comoving wavenumber kp associated with the "filtering scale" (note the 
change in convention from ir/kj to 2n/kF)- The wavenumber kp is related to 
the Jeans wavenumber kj by 



1 r , D{t') + 2H{t')D{t') r dt" 



klit) D 

where D{t) is the linear growth factor. At high redshift (where il^ 1), this 
relation simplifies to [153] 

' ^" -l-J'^ ) . (64) 




kUt) a Jo kj{a') 

Then the relationship between the linear overdensity of the dark matter Sdm 
and the linear overdensity of the baryons 5b, in the limit of small k, can be 
written as [150] 

^ = '-|-°('''- (»^> 

Linear theory specifies whether an initial perturbation, characterized by 
the parameters k, Sdm,i, <5b,i and ti, begins to grow. To determine the mini- 
mum mass of nonlinear baryonic objects resulting from the shell-crossing and 
virialization of the dark matter, we must use a different model which exam- 
ines the response of the gas to the gravitational potential of a virialized dark 
matter halo. 



3.6 Formation of Nonlinear Objects 

3.7 Spherical Collapse 

Let us consider a spherically symmetric density or velocity perturbation of the 
smooth cosmological background, and examine the dynamics of a test particle 
at a radius r relative to the center of symmetry. Birkhoff's (1923) [44] the- 
orem implies that we may ignore the mass outside this radius in computing 
the motion of our particle. We further find that the relativistic equations of 
motion describing the system reduce to the usual Friedmann equation for the 
evolution of the scale factor of a homogeneous Universe, but with a density 
parameter i? that now takes account of the additional mass or peculiar veloc- 
ity. In particular, despite the arbitrary density and velocity profiles given to 
the perturbation, only the total mass interior to the particle's radius and the 
peculiar velocity at the particle's radius contribute to the effective value of i?. 
We thus find a solution to the particle's motion wliic;h describes its departure 
from the background Hubble flow and its subsequent collapse or expansion. 
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This solution holds until our particle crosses paths with one from a different 

radius, which happens rather late for most initial profiles. 

As with the Friedmann equation for a smooth Universe, it is possible to 
reinterpret the problem into a Newtonian form. Here we work in an inertial 
(i.e. non-comoving) coordinate system and consider the force on the particle as 
that resulting from a point mass at the origin (ignoring the possible presence 
of a vacuum energy density): 

d'r GM 

where G is Newton's constant, r is the distance of the particle from the center 

of the spherical perturbation, and M is the total mass within that radius. As 
long as the radial shells do not cross each other, the mass M is constant in 
time. The initial density profile determines M, while the initial velocity profile 
determines dr/dt at the initial time. As is well-known, there are three branches 
of solutions: one in which the particle turns around and collapses, another in 
which it reaches an infinite radius with some asymptotically positive velocity, 
and a third intermediate case in which it reachs an infinite radius but with a 
velocity that approaches zero. These cases may be written as [164]: 



r = A{cosrj — 1) 
t = B{r] — sin rj) 



r = Arf/2 



} 



Closed (0 < 7/ < 27r) (67) 



Flat (0 < 7? < oo) (68) 



r 



t = Bif/Q 

where = GMB'^ apphes in all cases. All three solutions have = QGMt^ /2 
as t goes to zero, which matches the linear theory expectation that the per- 
turbation amplitude get smaller as one goes back in time. In the closed case, 
the shell turns around at time ttB and radius 2A and collapses to zero radius 
at time 2-kB. 

We are now faced with the problem of relating the spherical collapse pa- 
rameters A, B, and M to the linear theory density perturbation S [283]. We do 
this by returning to the equation of motion. Consider that at an early epoch 
(i.e. scale factor <C 1), we arc given a spherical patch of miiform overden- 
sity 6i (the so-called 'top-hat' perturbation). If ]7 is essentially unity at this 
time and if the perturbation is pure growing mode, then the initial velocity is 
radially inward with magnitude 5iH{ti)r /3, where H{ti) is the Hubble con- 
stant at the initial time and r is the radius from the center of the sphere. This 
can be easily seen from the continuity equation in spherical coordinates. The 
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equation of motion (in noncomoving coordinates) for a particle beginning at 
radius rj is simply 

d^r GM Ar 

where M = {4:n/3)rf pi{l+Si) and pi is the background density of the Universe 
at time ti. We next define the dimensionless radius x = rai/ri and rewrite 
equation (70) as 

Our initial conditions for the integration of this orbit are 



(71) 



x{ti) = ai 



Hoai I 1 - — 



'a 



(72) 
(73) 



where H{ti) = Ho[f2m/o,^{ti) + (1 — ^m)]^^"^ is the Hubble parameter for a 
flat Universe at a a cosmic time ti. Integrating equation (71) yields 



1 / dx 



iJ 2 V dt 



il + Si) + QAX^ + K, 



(74) 



where _ftr is a constant of integration. Evaluating this at the initial time and 
dropping terms of 0{ai) (but 6i ~ Ui, so we keep ratios of order unity), we 
flnd 

K = --^nm + flk. (75) 

If K is sufficiently negative, the particle will turn-around and the sphere will 
collapse at a time 



Z'""" o -1/2 

HotcoU = 2 / da {f2m/a + K + QaO^) , 
Jo 



(76) 



where amax is the value of a which sets the denominator of the integral to 
zero. 

For the case of yl = 0, wc can determine the spherical collapse parameters 
A and B. K > {K < 0) produces an open (closed) model. Comparing 
coefficients in the energy equations [eq. (74) and the integration of (66)], one 
finds 



A = 



2a,; 



55^ 
3a,; 



B 



2Ha 



-3/2 



(77) 



(78) 



where flk — ^ ^ ^m- In particular, in an i? = 1 Universe, where 1 + z ~ 
{3Hot/2)~'^/^, we find that a sheh collapses at redshift I + Zc = 0.5929(5j/ai, 
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or in other words a shell collapsing at redshift had a linear overdensity 
extrapolated to the present day of (5o 1.686(1 + z^). 

While this derivation has been for spheres of constant density, we may 
treat a general spherical density profile (5i(r) up until shell crossing [164]. A 
particular radial shell evolves according to the mass interior to it; therefore, 
we define the average overdensity 



so that we may use 8i in place of 5i in the above formulae. If (5i is not monoton- 
ically decreasing with i?, then the spherical top-hat evolution of two different 
radii will predict that they cross each other at sonic late time; this is known 
as shell crossing and signals the breakdown of the solution. Even well-behaved 
profiles will produce shell crossing if shells are allowed to collapse to r = 
and then reexpand, since these expanding shells will cross infalling shells. In 
such a case, first-time infalling shells will never be affected prior to their turn- 
around; the more complicated behavior after turn-around is a manifestation 
of virialization. While the end state for general initial conditions cannot be 
predicted, various results are known for a self-similar collapse, in which (5(r) 
is a power-law [132, 40], as well as for the case of secondary infall models 
[156, 165, 181]. 

3.8 Halo Properties 

The small density fluctuations evidenced in the CMB grow over time as de- 
scribed in the previous subsection, until the perturbation 5 becomes of order 
unity, and the full non-linear gravitational problem must be considered. The 
dynamical collapse of a dark matter halo can be solved analytically only in 
cases of particular symmetry. If we consider a region which is much smaller 
than the horizon cH~^, then the formation of a halo can be formulated as 
a problem in Newtonian gravity, in some cases with minor corrections com- 
ing from General Relativity. The simplest case is that of spherical symmetry, 
with an initial (f = <C io) top- hat of uniform overdensity Si inside a sphere 
of radius R. Although this model is restricted in its direct applicability, the 
results of spherical collapse have turned out to be surprisingly useful in un- 
derstanding the properties and distribution of halos in models based on cold 
dark matter. 

The collapse of a spherical top-hat perturbation is described by the New- 
tonian equation (with a correction for the cosmological constant) 




(79) 



GM 



(80) 



where r is the radius in a fixed (not comoving) coordinate frame, Hq is the 
present-day Hubble constant, M is the total mass enclosed within radius r. 
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and the initial velocity field is given by the Hubble flow dr/dt = H{t)r. The 
enclosed 5 grows initially as 5]^ = SiD(t)/ D(ti), in accordance with linear 
theory, but eventually S grows above S^.U the mass shell at radius r is bound 
(i.e., if its total Newtonian energy is negative) then it reaches a radius of max- 
imum expansion and subsequently collapses. As demonstrated in the previous 
section, at the moment when the top-hat collapses to a point, the overdensity 
predicted by linear theory is 5l = 1.686 in the Einstein-de Sitter model, with 
only a weak dependence on f2m and i7/i. Thus a tophat collapses at redshift 
z if its linear overdensity extrapolated to the present day (also termed the 
critical density of collapse) is 



^crit(^) = 



1.686 



(81) 



where we set D{z = 0) = 1. 

Even a slight violation of the exact symmetry of the initial perturbation 
can prevent the tophat from collapsing to a point. Instead, the halo reaches a 
state of virial equilibrium by violent relaxation (phase mixing). Using the virial 
theorem U = ~2K to relate the potential energy U to the kinetic energy K in 
the final state (implying that the virial radius is half the turnaround radius - 
where the kinetic energy vanishes), the final overdensity relative to the critical 
density at the collapse redshift is Ac = IStt^ ~ 178 in the Einstein-de Sitter 
model, modified in a Universe with + = 1 to the fitting formula (Bryan 
& Norman 1998 [71]) 

Z\c = IStt^ + 82d - 39d^ , (82) 
where d = — 1 is evaluated at the collapse redshift, so that 

A halo of mass M collapsing at redshift z thus has a virial radius 



0.784 



M 



10^ /i-i Mq 
and a corresponding circular velocity. 



^1/3 




Ac ] 


-1/3 






187r2 





10 



h-^ kpc , (84) 



Vc = 



GM 



1/2 



= 23.4 



M 



108 /i-i 



^1/3 




Ac 1 








187r2 





1/2 



km s 
(85) 

In these expressions we have assumed a present Hubble constant written in 
the form Hq = 100 /i km s~^Mpc~^. We may also define a virial temperature 



T., = ^ = 1.98x10^ 
2k 



(-) 
VO.6/ 



M 



108 /i-i Mp 



^2/3 




Ac ] 






Qz 


187r2 





(86) 



K, 
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where /x is the mean molecular weight and nip is the proton mass. Note that 
the value of /i depends on the ionization fraction of the gas; for a fully ionized 
primordial gas /U = 0.59, while a gas with ionized hydrogen but only singly- 
ionized helium has = 0.61. The binding energy of the halo is approximately^ 

^- = 2^ = ^-^^^^° ( ios/.-^mJ [nimrA (^j'^ 

(87) 

Note that the binding energy of the baryons is smaller by a factor equal to 
the baryon fraction Qh / fim ■ 

Although spherical collapse captures some of the physics governing the 
formation of halos, structure formation in cold dark matter models procedes 
hierarchically. At early times, most of the dark matter is in low-mass halos, 
and these halos continuously accrete and merge to form high-mass halos. Nu- 
merical simulations of hierarchical halo formation indicate a roughly universal 
spherically-averaged density profile for the resulting halos (Navarro, Prenk, & 
White 1997, hereafter NFW [266]), though with considerable scatter among 
different halos (e.g., [72]). The NFW profile has the form 

where x = r/rvii, and the characteristic density Sc is related to the concen- 
tration parameter cn by 

A „3 

(89) 



3 ln(l-hcN)-CN/(l + CN) ■ 



The concentration parameter itself depends on the halo mass M , at a given 
redshift z [377]. 

More recent N-body simulations indicate deviations from the original NFW 
profile; for details and refined fitting formula see [268]. 



4 Nonlinear Growth 

4.1 The Abundance of Dark Matter Halos 

In addition to characterizing the properties of individual halos, a critical pre- 
diction of any theory of structure formation is the abundance of halos, i.e. the 
number density of halos as a function of mass, at any redshift. This prediction 
is an important step toward inferring the abundances of galaxies and galaxy 
clusters. While the number density of halos can be measured for particular 

® Tlie coefHcient of 1/2 in equation (87) would be exact for a singular isothermal 
sphere, p(r) oc 1/r^. 
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cosmologies in numerical simulations, an analytic model helps us gain physical 
understanding and can bo used to explore the dependence of abundances on 
all the cosmological parameters. 

A simple analytic model which successfully matches most of the numerical 
simulations was developed by Press & Schcchtcr (1974) [291]. The model is 
based on the ideas of a Gaussian random field of density perturbations, linear 
gravitational growth, and spherical collapse. To determine the abundance of 
halos at a rcdshift z. wo use (5m, the density field smoothed on a mass scale 
M, as defined in §3.3. Since (5m is distributed as a Gaussian variable with 
zero mean and standard deviation cr(M) [which depends only on the present 
linear power spectrum, see equation (17)], the probability that 5m is greater 
than some (5 equals 



The fundamental ansatz is to identify this probability with the fraction of dark 
matter particles which are part of collapsed halos of mass greater than M, at 
redshift z. There are two additional ingredients: First, the value used for E is 
'5crit(-z) given in equation (81), which is the critical density of collapse found 
for a spherical top-hat (extrapolated to the present since fT(M) is calculated 
using the present power spectrum); and second, the fraction of dark matter 
in halos above M is multiplied by an additional factor of 2 in order to ensure 
that every particle ends up as part of some halo with M > 0. Thus, the final 
formula for the mass fraction in halos above M at redshift z is 



This ad-hoc factor of 2 is necessary, since otherwise only positive fluctu- 
ations of (5m would be included. Bond et al. (1991) [52] found an alternate 
derivation of this correction factor, using a different ansatz. In their derivation, 
the factor of 2 has a more satisfactory origin, namely the so-called "cloud-in- 
cloud" problem: For a given mass M, even if ^m is smaller than (5crit(^); it 
is possible that the corresponding region lies inside a region of some larger 
mass Mi > M, with > ^cvxtiz). In this case the original region should 
be counted as belonging to a halo of mass Mj,. Thus, the fraction of parti- 
cles which arc part of collapsed halos of mass greater than M is larger than 
the expression given in equation (90). Bond et al. showed that, under cer- 
tain assumptions, the additional contribution results precisely in a factor of 2 
correction. 

Differentiating the fraction of dark matter in halos above M yields the 
mass distribution. Letting dn be the comoving number density of halos of 
mass between M and M -\- dM, we have 




(90) 




(91) 




(92) 
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Fig. 13. Fraction of baryons that assembled into dark matter haJos with a virial 
temperature of Tvir > lO^'K as a function of redshift. These baryons are above the 
temperature threshold for gas cooling and fragmentation via atomic transitions. 

After reionization the temperature barrier for star formation in galaxies is raised 
because the photo-ionized intergalactic medium is already heated to ^ lO'^K and it 
can condense only into halos with Tvir > lO^K. 

where Vc = Scrii{z)/(T{M) is the number of standard deviations which the 
critical collapse overdensity represents on mass scale M. Thus, the abundance 
of halos depends on the two functions a{M) and Sciit{z), each of which depends 
on the energy content of the Universe and the values of the other cosmological 
parameters. Since recent observations confine the standard set of parameters 
to a relatively narrow range, we illustrate the abundance of halos and other 
results for a single set of parameters: = 0.3, Ha = 0.7, fib = 0.045, 
as = 0.9, a primordial power spectrum index n = 1 and a Hubble constant 
h = 0.7. 

Figure 14 shows (t{M) and Scrit{z), with the input power spectrum com- 
puted from Eisenstein & Hu (1999) [118]. The sohd line is a{M) for the cold 
dark matter model with the parameters specified above. The horizontal dot- 
ted lines show the value of (5crit(z) at z = 0, 2, 5, 10, 20 and 30, as indicated in 
the figure. From the intersection of these horizontal lines with the solid line 
we infer, e.g., that at2; = 5al — a fluctuation on a mass scale of 2 x 1O''M0 
will collapse. On the other hand, at z = 5 collapsing halos require a 2 — cr fluc- 
tuation on a mass scale of 3 x 10^^ Mq, since a{M) on this mass scale equals 
about half of (5crit(^ = 5). Since at each redshift a fixed fraction (31.7%) of the 
total dark matter mass lies in halos above the 1 — cr mass. Figure 14 shows 
that most of the mass is in small halos at high redshift, but it continuously 
shifts toward higher characteristic halo masses at lower redshift. Note also 
that cr(M) flattens at low masses because of the changing shape of the power 
spectrum. Since cr ^ oo as M ^ 0, in the cold dark matter model all the 
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dark matter is tied up in halos at all redshifts, if sufficiently low-mass halos 
are considered. 
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Fig. 14. Mass fluctuations and collapse thresholds in cold dark matter models. 
The horizontal dotted lines show the value of the extrapolated collapse overdensity 
5crit(z) at the indicated redshifts. Also shown is the value of a{M) for the cos- 
mological parameters given in the text (solid curve), as well as a{M) for a power 
spectrum with a cutoff below a mass M = 1.7 x 10* Mq (short-dashed curve), or 
M = 1.7 X 10"Mq (long-dashed curve). The intersection of the horizontal lines 
with the other curves indicate, at each redshift z, the mass scale (for each model) 
at which a 1 — a fluctuation is just collapsing at z (see the discussion in the text). 

Also shown in Figure 14 is the effect of cutting off the power spectrum on 
small scales. The short-dashed curve corresponds to the case where the power 
spectrum is set to zero above a comoving wavenumber k = lOMpc""'^, which 
corresponds to a mass M = 1.7 x IO^Mq. The long-dashed curve corresponds 
to a more radical cutoff above k = lMpc~^, or below M = 1.7 x lO^^Af©. A 
cutoff severely reduces the abundance of low-mass halos, and the finite value 
of a{M — 0) implies that at all redshifts some fraction of the dark matter 
does not fall into halos. At high redshifts where 5crit{z) ^ a{M — 0), all 
halos are rare and only a small fraction of the dark matter lies in halos. In 
particular, this can affect the abundance of halos at the time of reionization, 
and thus the observed limits on reionization constrain scenarios which include 
a small-scale cutoff in the power spectrum [21]. 

In figures 15 - 18 we show explicitly the properties of collapsing halos 
which represent 1 — cr, 2 — cr, and 3 — cr fluctuations (corresponding in all cases 
to the curves in order from bottom to top), as a function of redshift. No cutoff 
is applied to the power spectrum. Figure 15 shows the halo mass. Figure 16 
the virial radius. Figure 17 the virial temperature (with /i in equation (86) 
set equal to 0.6, although low temperature halos contain neutral gas) as well 
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as circular velocity, and Figure 18 shows the total binding energy of these 
halos. In figures 15 and 17, the dotted curves indicate the minimum virial 
temperature required for efficient cooling with primordial atomic species only 
(upper curve) or with the addition of molecular hydrogen (lower curve). Figure 
18 shows the binding energy of dark matter halos. The binding energy of the 
baryons is a factor ~ Qb/ ~ 15% smaller, if they follow the dark matter. 
Except for this constant factor, the figure shows the minimum amount of 
energy that needs to be deposited into the gas in order to unbind it from 
the potential well of the dark matter. For example, the hydrodynamic energy 
released by a single supernovae, ~ 10^"'^ erg, is sufficient to unbind the gas in 
all 1 — fj halos at z > 5 and in all 2 — cr halos at z > 12. 
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Fig. 15. Characteristic properties of collapsing halos: Halo mass. The solid curves 
show the mass of collapsing halos which correspond to 1 — a, 2 — cr, and 3 — cr fluctua- 
tions (in order from bottom to top). The dotted curves show the mass corresponding 
to the minimum temperature required for efficient cooling witli primordial atomic 
species only (upper curve) or with the addition of molecular hydrogen (lower curve). 

At z = 5, the halo masses which correspond to 1 — cr, 2 — cr, and 3 — cr 
fluctuations are 1.8 x IO'^Mq, 3.0 x IQ^^Mq, and 7.0 x IO^Mq, respectively. 
The corresponding virial temperatures are 2.0 x lO'^ K, 2.8 x 10^ K, and 
2.3 X 10^ K. The equivalent circular velocities are 7.5 km s^^, 88 km s^^, and 
250 km s^^. At z = 10, the 1 — cr, 2 — cr, and 3 — cr fluctuations correspond 
to halo masses of 1.3 X IO^Mq, 5.7 X IO'^Mq, and 4.8 x lO^M©, respectively. 
The corresponding virial temperatures are 6.2 K, 7.9 x 10^ K, and 1.5 x 
10^ K. The equivalent circular velocities are 0.41 km s^^, 15 km s~^, and 65 
km s^^. Atomic cooling is efficient at Tvir > 10^ K, or a circular velocity 
Vc > 17 km s^^. This corresponds to a 1.2 — cr fluctuation and a halo mass of 
2.1 x IO^Mq at z = 5, and a 2.1 — cr fluctuation and a halo mass of 8.3 x 10 ''M© 
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Fig. 16. Characteristic properties of collapsing halos: Halo virial radius. The curves 
show the virial radius of collapsing halos which correspond to 1 — cr, 2 — cr, and 3 — u 
fluctuations (in order from bottom to top). 
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Fig. 17. Characteristic properties of collapsing halos: Halo virial temperature and 
circular velocity. The solid curves show the virial temperature (or, equivalently, the 
circular velocity) of collapsing halos which correspond to 1 — a, 2 — a, and 3 — a 
fluctuations (in order from bottom to top). The dotted curves show the minimum 
temperature required for efficient cooling with primordial atomic species only (upper 
curve) or with the addition of molecular hydrogen (lower curve). 
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Fig. 18. Characteristic properties of collapsing halos: Halo binding energy. The 
curves show the total binding energy of collapsing halos which correspond to 1 — a, 
2 — cr, and 3 — a fluctuations (in order from bottom to top). 



at z = 10. Molecular hydrogen provides efficient cooling down to Tvir ~ 300 
K, or a circular velocity Vc ~ 2.9 km s~^. This corresponds to a 0.81 — a 
fluctuation and a halo mass of 1.1 x lO^M© at z = 5, and a 1.4 — cr fluctuation 
and a halo mass of 4.3 x lO^M© at z = 10. 

In Figure 19 we show the halo mass function dn/d\n(M) at several different 
redshifts: z — (solid curve), z — 5 (dotted curve), z = 10 (short-dashed 
curve) , z = 20 (long-dashed curve) , and z = 30 (dot-dashed curve) . Note that 
the mass function does not decrease monotonically with redshift at all masses. 
At the lowest masses, the abundance of halos is higher at z > than at z = 0. 

4.2 The Excursion-Set (Extended Press-Schechter) Formalism 

The usual Press-Schechter formalism makes no attempt to deal with the cor- 
relations between halos or between different mass scales. In particular, this 
means that while it can generate a distribution of halos at two different epochs, 
it says nothing about how particular halos in one epoch are related to those 
in the second. We therefore would like some method to predict, at least sta- 
tistically, the growth of individual halos via accretion and mergers. Even re- 
stricting ourselves to spherical collapse, such a model must utilize the full 
spherically-averaged density profile around a particular point. The potential 
correlations between the mean overdensities at different radii make the sta- 
tistical description substantially more difficult. 

The excursion set formalism (Bond et al. 1991 [52]) seeks to describe the 
statistics of halos by considering the statistical properties of S{R), the average 
overdensity within some spherical window of characteristic radius R, as a 
function of R. While the Press-Schechter model depends only on the Gaussian 
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Fig. 19. Halo mass function at several redshifts: z = (solid curve), z = 5 (dotted 
curve), z = 10 (short-dashed curve), z = 20 (long-dashed curve), and 2 = 30 (dot- 
dashed curve). 

distribution of i5 for one particular R, the excursion set considers all R. Again 
the connection between a value of the linear regime d and the final state is 
made via the spherical collapse solution, so that there is a critical value Sc{z) 
of d which is required for collapse at a redshift z. 

For most choices of window function, the functions 6{R) are correlated 
from one R to another such that it is prohibitively difficult to calculate the 
desired statistics directly [although Monte Carlo realizations are possible [52]]. 
However, for one particular choice of a window function, the correlations be- 
tween different R greatly simplify and many interesting quantities may be 
calculated [52, 212]. The key is to use a /c-space top-hat window function, 
namely Wk = 1 for all k less than some critical kc and Wk = ioi k > kc- 
This filter has a spatial form of W{r) cx ji{kcr)/kcr, which implies a volume 
6n'^/k^ or mass Qn'^pb/k^. The characteristic radius of the filter is ~ k~^ , as 
expected. Note that in real space, this window function converges very slowly, 
due only to a sinusoidal oscillation, so the region under study is rather poorly 
localized. 

The great advantage of the sharp fc-space filter is that the difference at 
a given point between 5 on one mass scale and that on another mass scale 
is statistically independent from the value on the larger mass scale. With a 
Gaussian random field, each Sk is Gaussian distributed independently from 
the others. For this filter. 



meaning that the overdensity on a particular scale is simply the sum of the 
random variables Sk interior to the chosen kc- Consequently, the difference 




(93) 
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between the S{M) on two mass scales is just the sum of the Sk in the spher- 
ical shell between the two fee, which is independent from the sum of the S). 
interior to the smaller kc- Meanwhile, the distribution of d{M) given no prior 
information is still a Gaussian of mean zero and variance 

a\M) = ^[ dkk^P{k). (94) 

^T"" Jk<kc{M) 

If we now consider 5 as a function of scale fee, we see that we begin from 
^ = at fcc = (M = oo) and then add independently random pieces as 
kc increases. This generates a random walk, albeit one whose stepsize varies 
with kc- We then assume that, at redshift z, a given function 5{kc) represents 
a collapsed mass M corresponding to the kc where the function first crosses 
the critical value 5c{z). With this assumption, we may use the properties of 
random walks to calculate the evolution of the mass as a function of redshift. 

It is now easy to rederive the Press-Schcchtcr mass function, including the 
previously unexplained factor of 2 [52, 212, 382]. The fraction of mass elements 
included in halos of mass less than M is just the probability that a random 
walk remains below 5c{z) for all kc less than Kc, the filter cutoff appropriate 
to M . This probability must be the complement of the sum of the probabilities 
that: (a) 5{Kc) > Sc{z); or that (b) S{Kc) < 6c{z) but S{k'J > 5c{z) for some 
k'^ < Kc- But these two cases in fact have equal probability; any random walk 
belonging to class (a) may be reflected around its first upcrossing of 6c{z) to 
produce a walk of class (b), and vice versa. Since the distribution of 6{Kc) is 
simply Gaussian with variance (T^(Af), the fraction of random walks falling 
into class (a) is simply {1 / V2na^) f°^f^y^ dS cxp{—S^/2a^{M)}. Hence, the 
fraction of mass elements included in halos of mass less than M at redshift z 
is simply 

1 f°° 

F{< M) = 1 - 2 X —= / dS exp{-(5V2cr^M)} (95) 

V 2770-2 JSc{z) 

which may be diff'erentiated to yield the Press-Schechter mass function. We 

may now go further and consider how halos at one redshift arc related to those 
at another redshift. If we are given that a halo of mass M2 exists at redshift 
Z2, then we know that the random function 5{kc) for each mass element within 
the halo first crosses S{z2) at kc2 corresponding to M2- Given this constraint, 
we may study the distribution of kc where the function 6{kc) crosses other 
thresholds. It is particularly easy to construct the probability distribution for 
when trajectories first cross some Sc{zi) > Sc{z2) (implying zi > Z2); clearly 
this occurs at some kci > kc2 ■ This problem reduces to the previous one if we 
translate the origin of the random walks from {kc, 5) = (0, 0) to {kc2, 5c{z2))- 
We therefore find the distribution of halo masses Mi that a mass element 
finds itself in at redshift zi given that it is part of a larger halo of mass M2 
at a later redshift Z2 is [52, 55]) 
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dP 



(Mi,Zl|M2,Z2) 



(96) 




2 5c{zi) - 5^{z2) da {Ml) 

7r[tT2(Mi)-(T2(M2)]3/2 dMi 



exp 



2[a2(Mi)-a2(M2)] 



[<5c(^l)- -5,(22)]' 



} 



This may be rewritten as saying that the quantity 



V = 



Sc{zi) - Sc{z2) 



(97) 



is distributed as the positive half of a Gaussian with unit variance; equation 
(97) may be inverted to find Mi{v). 

We seek to interpret the statistics of these random walks as those of merg- 
ing and accreting halos. For a single halo, we may imagine that as we look back 
in time, the object breaks into ever smaller pieces, similar to the branching of 
a tree. Equation (96) is the distribution of the sizes of these branches at some 
given earlier time. However, using this description of the ensemble distribution 
to generate Monte Carlo realizations of single merger trees has proven to be 
difficult. In all cases, one recursively steps back in time, at each step breaking 
the final object into two or more pieces. An elaborate scheme (Kauffmann 
& White 1993 [195]) picks a large number of progenitors from the ensemble 
distribution and then randomly groups them into sets with the correct total 
mass. This generates many (hundreds) possible branching schemes of equal 
likelihood. A simpler scheme (Lacey & Cole 1993 [212]) assumes that at each 
time step, the object breaks into two pieces. One value from the distribution 
(96) then determines the mass ratio of the two branchs. 

One may also use the distribution of the ensemble to derive some additional 
analytic results. A useful example is the distribution of the epoch at which 
an object that has mass M2 at redshift Z2 has accumulated half of its mass 
[212]. The probability that the formation time is earlier than Zi is equal to the 
probability that at redshift zi a progenitor whose mass exceeds M2/2 exists: 



where dP/dM is given in equation (96). The factor of M2/M corrects the 
counting from mass weighted to number weighted; each halo of mass M2 can 
have only one progenitor of mass greater than M2/2. Differentiating equation 
(98) with respect to time gives the distribution of formation times. This an- 
alytic form is an excellent match to scale-free N-body simulations [213]. On 
the other hand, simple Monte Carlo implementations of equation (96) produce 
formation redshifts about 40% higher [212]. As there may be correlations be- 
tween the various branches, there is no unique Monte Carlo scheme. 

Numerical tests of the excursion set formalism are quite encouraging. Its 
predictions for merger rates arc in very good agreement with those measured 
in scale-free N-body simulations for mass scales down to around 10% of the 




(98) 
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nonUnear mass scale (that scale at which <jm = 1 ) ) and distributions of for- 
mation times closely match the analytic predictions [213]. The model appears 
to be a promising method for tracking the merging of halos, with many ap- 
plications to cluster and galaxy formation modeling. In particular, one may 
use the formalism as the foundation of semi-analytic galaxy formation models 
[196]. The excursion set formalism may also be used to derive the correlations 
of halos in the nonlinear regime [258] . 



4.3 Response of Baryons to Nonlinear Dark Matter Potenticds 

The dark matter is assumed to be cold and to dominate gravity, and so its 
collapse and virialization proceeds unimpeded by pressure effects. In order to 
estimate the minimum mass of baryonic objects, we must go beyond linear 
perturbation theory and examine the baryonic mass that can accrete into the 
final gravitational potential well of the dark matter. 

For this purpose, we assume that the dark matter had already virialized 
and produced a gravitational potential 0(r) at a redshift z^ir (with </> ^ at 
large distances, and (/) < inside the object) and calculate the resulting over- 
density in the gas distribution, ignoring cooling (an assumption justified by 
spherical collapse simulations which indicate that cooling becomes important 
only after virialization; see Haiman et al. 1996 [168]). 

After the gas settles into the dark matter potential well, it satisfies the 
hydrostatic equilibrium equation, 

Vpb = -PbV</> (99) 

where p\, and pb are the pressure and mass density of the gas. At z < 100 
the gas temperature is decoupled from the CMB, and its pressure evolves 
adiabatically (ignoring atomic or molecular cooling). 



Ph \PhJ 



5/3 

(100) 



where a bar denotes the background conditions. We substitute equation (100) 
into (99) and get the solution, 



Ph 



5 kT J 



where T = phl^rnp/{kph) is the backgromid gas temperature. If we define 
Tvir = —^mp(j)/k as the virial temperature for a potential depth —cj), then the 
overdensity of the baryons at the virialization redshift is 
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This solution is approximate for two reasons: (i) we assumed that the gas is 
stationary throughout the entire region and ignored the transitions to infall 
and the Hubble expansion at the interface between the collapsed object and 
the background intergalactic medium (henceforth IGM), and (ii) we ignored 
entropy production at the virialization shock surrounding the object. Never- 
theless, the result should provide a better estimate for the minimum mass of 
collapsed baryonic objects than the Jeans mass does, since it incorporates the 
nonlinear potential of the dark matter. 

We may define the threshold for the collapse of baryons by the criterion 
that their mean overdensity, 5b, exceeds a value of 100, amounting to > 50% 
of the baryons that would assemble in the absence of gas pressure, according 
to the spherical top-hat collapse model. Equation (102) then implies that 
Tvir > 17.2 T. 

As mentioned before, the gas temperature evolves at z < 160 according to 
the relation f 170[(l + 2)/100]2 K. This implies that baryons are overdense 
by > 100 only inside halos with a virial temperature Tvir > 2.9 x 10^ [(1 -|- 
z)/100]^ K. Based on the top-hat model, this implies a minimum halo mass 
for baryonic objects of 

M.. = .0..0»(^)-"^(i|S)-"'(l±£)"'M., a03, 

where we consider sufficiently high redshifts so that ~ 1. This minimum 
mass is coincident ally almost identical to the naive Jeans mass calculation 
of linear theory in equation (62) despite the fact that it incorporates shell 
crossing by the dark matter, which is not accounted for by linear theory. Unlike 
the Jeans mass, the minimum mass depends on the choice for an overdensity 
threshold [taken arbitrarily as > 100 in equation (103)]. To estimate the 
minimum halo mass which produces any significant accretion we set, e.g., 
Sh = 5, and get a mass which is lower than Mmin by a factor of 27. 

Of course, once the first stars and quasars form they heat the surrounding 
IGM by either outflows or radiation. As a result, the Jeans mass which is rele- 
vant for the formation of new objects changes [148, 152]). The most dramatic 
change occurs when the IGM is photo-ionized and is consequently heated to 
a temperature of ~ (1-2) x 10'' K. 

5 Fragmentation of the First Gaseous Objects to Stars 

5.1 Star Formation 

As mentioned in the preface, the fragmentation of the first gaseous objects is 
a well-posed physics problem with well specified initial conditions, for a given 
power-spectrum of primordial density fluctuations. This problem is ideally 
suitc;d for three-dimensional computer simulations, since it cannot be reliably 
addressed in idealized ID or 2D geometries. 
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Recently, two groups have attempted detailed 3D simulations of the forma- 
tion process of the first stars in a halo of ~ lO^AfQ by following the dynamics 
of both the dark matter and the gas components, including H2 chemistry and 
cooling. Bromm, Coppi, & Larson (1999) [57] have used a Smooth Particle 
Hydrodynamics (SPH) code to simulate the collapse of a top-hat overden- 
sity with a prescribed solid-body rotation (corresponding to a spin parameter 
A = 5%) and additional small perturbations with P{k) cx k^^ added to the 
top-hat profile. Abel et al. (2002) [5] isolated a high-density filament out of a 
larger simulated cosmological volume and followed the evolution of its density 
maximum with exceedingly high resolution using an Adaptive Mesh Refine- 
ment (AMR) algorithm. 




T [K] 



Fig. 20. Cooling rates as a function of temperature for a primordial gas composed 
of atomic hydrogen and helium, as well as molecular hydrogen, in the absence of 
any external radiation. We assume a hydrogen number density uh ~ 0.045 cm~^, 
corresponding to the mean density of virialized halos at 2: = 10. The plotted quantity 
A/n^ is roughly independent of density (unless uh > 10 cm~^), where A is the 
volume cooling rate (in erg/sec/cm^). The solid line shows the cooling curve for 
an atomic gas, with the characteristic peaks due to coUisional excitation of HI 
and He2. The dashed line shows the additional contribution of molecular cooling, 
assuming a molecular abundance equal to 1% of uh- 

The generic results of Bromm et al. (1999 [57]; see also Bromm 2000 [58]) 
are illustrated in Figure 21. The collapsing region forms a disk which fragments 
into many clumps. The clumps have a typical mass ~ 10^-10"^ Mq. This mass 
scale corresponds to the Jeans mass for a temperature of ~ 500K and the 
density ~ 10'* cm~^ where the gas lingers because its cooling time is longer 
than its collapse time at that point (see Fig. 22). Each clump accretes mass 
slowly until it exceeds the Jeans mass and collapses at a roughly constant 
temperature (isothermally) due to H2 cooling that brings the gas to a fixed 
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temperature floor. The clump formation efficiency is high in this simulation 
due to the synchronized collapse of the overall top-hat perturbation. 




Fig. 21. Numerical results from Bromm et al. (1999) [57], showing gas properties 
at z — 31.2 for a eoUapsing slightly inhomogeneous top-hat region with a prescribed 
solid-body rotation, (a) Free electron fraction (by number) vs. hydrogen number 
density (in cm~^). At densities exceeding n ~ 10^ cm~^, recombination is very 
efficient, and the gas becomes almost completely neutral, (b) Molecular hydrogen 
fraction vs. number density. After a quick initial rise, the H2 fraction approaches 
the asymptotic value of / ~ 10~^, due to the H~ channel, (c) Gas temperature vs. 
number density. At densities below ~ 1 cm~^, the gas temperature rises because of 
adiabatic compression until it reaches the virial value of T^ir — 5000 K. At higher 
densities, cooling due to H2 drives the temperature down again, until the gas settles 
into a quasi-hydrostatic state at T ~ 500 K and n ~ 10** cm~^. Upon further 
compression due to accretion and the onset of gravitational collapse, the gas shows 
a further modest rise in temperature, (d) Jeans mass (in Mq) vs. number density. 
The Jeans mass reaches a value of Mj ~ 10"' Mq for the quasi- hydrostatic gas in 
the center of the potential well, and reaches the resolution limit of the simulation, 
Mres ^ 2OOM0, for densities close to n = 10® cm"^. 



Bromm (2000) [58] has simulated the collapse of one of the above- 
mentioned clumps with ~ IOOOMq and demonstrated that it does not tend 
to fragment into sub-components. Rather, the clump core of ~ IOOMq free- 
falls towards the center leaving an extended envelope behind with a roughly 
isothermal density proffie. At very high gas densities, thrcc-body reactions 
become important in the chemistry of H2. Omukai & Nishi (1998) [274] have 
included these reactions as well as radiative transfer and followed the collapse 
in spherical symmetry up to stellar densities. Radiation pressure from nuclear 
burning at the center is unlikely to reverse the infall as the stellar mass builds 
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Fig. 22. Gas and clump morphology z = 28.9 in the simulation of Bromm 
et al. (1999) [57]. Top row: The remaining gas in the diffuse phase. Bottom row: 
Distribution of clumps. The numbers next to the dots denote clump mass in units of 

Mq . Left panels: Facc-on view. Right panels: Edgc-on view. The length of the box is 
30 pc. The gas has settled into a flattened configuration with two dominant clumps 
of mass close to 20, DOOM©. During the subsequent evolution, the clumps survive 
without merging, and grow in mass only slightly by accretion of surrounding gas. 



up. These calculations indicate that each clump may end as a single mas- 
sive star; however, it is conceivable that angular momentum may eventually 
halt the collapsing cloud and lead to the formation of a binary stellar system 
instead. 

The Jeans mass, which is defined based on small fluctuations in a back- 
ground of uniform density, does not strictly apply in the context of collapsing 
gas cores. We can instead use a slightly modified critical mass known as the 
Bonnor-Ebert mass [53, 114]. For baryons in a background of uniform den- 
sity pb, perturbations are unstable to gravitational collapse in a region more 
massive than the Jeans mass. Instead of a uniform background, we consider 
a spherical, non-singular, isothermal, self-gravitating gas in hydrostatic equi- 
librium, i.e., a centrally-concentrated object which more closely resembles the 
gas cores found in the above-mentioned simulations. In this case, small fluc- 
tuations are unstable and lead to collapse if the sphere is more massive than 
the Bonnor-Ebert mass Mbe, given by the same expression the Jeans Mass 
but with a different coefficient (1.2 instead of 2.9) and with pf, denoting in 
this case the gas (volume) density at the surface of the sphere, 

1 / hT \ 

Mbe = 1.2— • (104) 

^ \GfimpJ 

In their simulation, Abel et al. (2000) [4] adopted the actual cosmological 
density perturbations as initial conditions. The simulation focused on the den- 
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sity peak of a filament within the IGM, and evolved it to very high densities 
(Fig. 23). Following the initial collapse of the filament, a clump core formed 
with ~ 200Mq, amounting to only ~ 1% of the virialized mass. Subsequently 
due to slow cooling, the clump collapsed subsonically in a state close to hy- 
drostatic equilibrium (see Fig. 24). Unlike the idealized top-hat simulation 
of Bromm et al. (2001) [59], the collapse of the different clumps within the 
filament is not synchronized. Once the first star forms at the center of the 
first collapsing clump, it is likely to affect the formation of other stars in its 
vicinity. 

As soon as nuclear burning sets in the core of the proto-star, the radiation 
emitted by the star starts to affect the infall of the surrounding gas towards 
it. The radiative feedback involves photo-dissociation of H2, Lya radiation 
pressure, and photo-evaporation of the accretion disk. Tan & McKee [357] 
studied these effects by extrapolating analytically the infall of gas from the 
final snapshot of the above resolution-limited simulations to the scale of a 
proto-star; they concluded that nuclear burning (and hence the feedback) 
starts when the proton-star accretes ~ SOM© and accretion is likely to be 
terminated when the star reaches ^ 200Mq. 




Fig. 23. Zooming in on the core of a star forming region with the Adaptive Mesh 
Refinement simulation of Abel et al. (2000) [4]. The panels show different length 
scales, decreasing clockwise by an order of magnitude between adjacent panels. Note 
the large dynamic range of scales which are being resolved, from 6 kpc (top left panel) 
down to 10,000 AU (bottom left panel). 



If the clumps in the above simulations end up forming individual very 
massive stars, then these stars will likely radiate copious amounts of ionizing 
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Fig. 24. Gas profiles from the simulation of Abel et al. (2000) [4]. The cell size on 
the finest grid corresponds to 0.024 pc, while the simulation box size corresponds to 
6.4 kpc. Shown are spherically- averaged mass-weighted profiles around the baryon 
density peak shortly before a well defined fragment forms (z — 19.1). Panel (a) 
shows the baryonic number density, enclosed gas mass in solar mass, and the local 
Bonnor-Ebert mass Mbe (see text). Panel (b) plots the molecular hydrogen fraction 
(by number) and the free electron fraction x. The H2 cooling time, tH2, the 
time it takes a sound wave to travel to the center, fcross, and the free-fall time 
tff = [37r/(32Gp)]^''^ arc given in panel (c). Panel (d) gives the temperature in K 
as a function of radius. The bottom panel gives the local sound speed, Cs (solid line 
with circles) , the rms radial velocities of the dark matter (dashed line) and the gas 
(dashed line with asterisks) as well as the rms gas velocity (solid line with square 
symbols). The vertical dotted line indicates the radius (~ 5 pc) at which the gas 
has reached its minimum temperature allowed by H2 cooling. The virial radius of 
the 5.6 X lO'^Mo halo is 106 pc. 

radiation [50, 370, 59] and expel strong winds. Hence, the stars will have a large 
effect on their interstellar environment, and feedback is likely to control the 
overall star formation efficiency. This efSciency is likely to be small in galactic 
potential wells which have a virial temperature lower than the temperature 
of photoionized gas, ~ lO^K. In such potential wells, the gas may go through 
only a single generation of star formation, leading to a "suicidal" population 
of massive stars. 

The final state in the evolution of these stars is uncertain; but if their mass 
loss is not too extensive, then they are likely to end up as black holes [50, 137]. 
The remnants may provide the seeds of quasar black holes [215]. Some of the 
massive stars may end their lives by producing gamma-ray bursts. If so then 
the broad-band afterglows of these bursts could provide a powerful tool for 
probing the epoch of reionization [214, 94]). There is no better way to end the 
dark ages than with 7-ray burst fireworks. 
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Where are the first stars or their remnants located today? The very first 
stars formed in rare high-cr peaks and hence are hkely to populate the cores 
of present-day galaxies [380]. However, the bulk of the stars which formed 
in low-mass systems at later times are expected to behave similarly to the 
coUisionless dark matter particles and populate galaxy halos [221]. 



5.2 The Mass Function of Stars 



Currently, we do not have direct observational constraints on how the first 
stars, the so-called Population III stars, formed at the end of the cosmic dark 
ages. It is, therefore, instructive to briefly summarize what we have learned 
about star formation in the present-day Universe, where theoretical reasoning 
is guided by a wealth of observational data (see [293] for a recent review). 

Population I stars form out of cold, dense molecular gas that is structured 
in a complex, highly inhomogeneous way. The molecular clouds are supported 
against gravity by turbulent velocity fields and pervaded on large scales by 
magnetic fields. Stars tend to form in clusters, ranging from a few hundred up 
to ~ 10^ stars. It appears likely that the clustered nature of star formation 
leads to complicated dynamics and tidal interactions that transport angular 
momentum, thus allowing the collapsing gas to overcome the classical centrifu- 
gal barrier [216]. The initial mass function (IMF) of Pop I stars is observed 
to have the approximate Salpeter form (e.g., [208]) 

dN 

oc , (105) 



dlogM 
where 

(106) 



-1.35 for M > 0.5Mq 

0.0 for 0.007 <M< 0.5Mq 



The lower cutoff in mass corresponds roughly to the opacity limit for frag- 
mentation. This limit reflects the minimum fragment mass, set when the rate 
at which gravitational energy is released during the collapse exceeds the rate 
at which the gas can cool (e.g., [298]). The most important feature of the 
observed IMF is that ~ IMq is the characteristic mass scale of Pop I star for- 
mation, in the sense that most of the mass goes into stars with masses close 
to this value. In Figure 25, we show the result from a recent hydro dynamical 
simulation of the collapse and fragmentation of a molecular cloud core [31, 32]. 
This simulation illustrates the highly dynamic and chaotic nature of the star 
formation process^. 

Th{^ nicital-rich chemistry, magnctohydrodynamics, and radiative transfer 
involved in present-day star formation is complex, and we still lack a com- 
prehensive theoretical framework that predicts the IMF from first principles. 
Star formation in the high redshift Universe, on the other hand, poses a the- 
oretically more tractable problem due to a number of simplifying features, 



^ See http:// www.ukaff.ac.uk/starcluster for an animation. 
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Fig. 25. A hydrodynamic simulation of the coUapse and fragmentation of a turbu- 
lent molecular cloud in the present-day Universe (from [32]). The cloud has a mass 
of 50Mq. The panels show the column density through the cloud, and span a scale 
of 0.4 pc across. Left: The initial phase of the collapse. The turbulence organizes 
the gas into a network of filaments, and decays thereafter through shocks. Right: 
A snapshot taken near the end of the simulation, after 1.4 initial free-fall times of 
2 X lO^yr. Fragmentation has resulted in ~ 50 stars and brown dwarfs. The star 
formation efficiency is ~ 10% on the scale of the overall cloud, but can be much 
larger in the dense sub-condensations. This result is in good agreement with what 
is observed in local star-forming regions. 

such as: (i) the initial absence of heavy metals and therefore of dust; and 
(ii) the absence of dynamically-significant magnetic fields, in the pristine gas 
left over from the big bang. The cooling of the primordial gas does then only 
depend on hydrogen in its atomic and molecular form. Whereas in the present- 
day interstellar medium, the initial state of the star forming cloud is poorly 
constrained, the corresponding initial conditions for primordial star forma- 
tion are simple, given by the popular ylCDM model of cosmological structure 
formation. We now turn to a discussion of this theoretically attractive and 
important problem. 

How did the first stars form? A complete answer to this question would 
entail a theoretical prediction for the Population III IMF, which is rather 
challenging. Let us start by addressing the simpler problem of estimating the 
characteristic mass scale of the first stars. As mentioned before, this mass 
scale is observed to be ~ ^Mq in the present-day Universe. 

Bromm & Loeb (2004) [67] carried out idealized simulations of the proto- 
stellar accretion problem and estimated the final mass of a Population III star. 
Using the smoothed particle hydrodynamics (SPH) method, they included the 
chemistry and cooling physics relevant for the evolution of metal-free gas (see 
[62] for details). Improving on earlier work [57, 62] by initializing the simula- 
tions according to the ylCDM model, they focused on an isolated overdense 
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region that corresponds to a 3cr— peak [67]: a halo containing a total mass of 
IO^Mq, and collapsing at a redshift Zvir — 20. In these runs, one high-density 
clump has formed at the center of the minihalo, possessing a gas mass of a 
few hundred solar masses. Soon after its formation, the clump becomes grav- 
itationally unstable and undergoes runaway collapse. Once the gas clump has 
exceeded a threshold density of 10^ cm~'^, it is replaced by a sink particle 
which is a collisionless point-like particle that is inserted into the simulation. 
This choice for the density threshold ensures that the local Jeans mass is 
resolved throughout the simulation. The clump (i.e., sink particle) has an ini- 
tial mass of Mci — 200Mq, and grows subsequently by ongoing accretion of 
surrounding gas. High-density clumps with such masses result from the chem- 
istry and cooling rate of molecular hydrogen, H2, which imprint characteristic 
values of temperature, T ~ 200 K, and density, n ^ 10'' cm""^, into the metal- 
free gas [62] . Evaluating the Jeans mass for these characteristic values results 
in Mj ~ a few x IQ^Mq, which is close to the initial clump masses found 
in the simulations. 




Fig. 26. Collapse and fragmentation of a primordial cloud (from [67]). Shown is 
the projected gas density at a redshift z ~ 21.5, briefly after gravitational runaway 
collapse has commenced in the center of the cloud. Left: The coarse-grained mor- 
phology in a box with linear physical size of 23.5 pc. At this time in the unrefined 
simulation, a high-density clump (sink particle) has formed with an initial mass of 
IG^ Mq. Right: The refined morphology in a box with linear physical size of 0.5 pc. 
The central density peak, vigorously gaining mass by accretion, is accompanied by 
a secondary clump. 

The high-density clumps are clearly not stars yet. To probe the subsequent 
fate of a clump, Bromm & Loeb (2004) [67] have re-simulated the evolution 
of the central clump with sufficient resolution to follow the collapse to higher 
densities. Figure 26 [right panel) shows the gas density on a scale of 0.5 pc. 
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Fig. 27. Accretion onto a primordial protostar (from [67]). The morphology of this 
accretion flow is shown in Fig. 26. Left: Accretion rate (in M© yr"'^) vs. time (in yr) 

since molecular core formation. Right: Mass of the central core (in Mq) vs. time. 
Solid line: Accretion history approximated as: M» oc t" ''^. Using this analytical 
approximation, we extrapolate that the protostellar mass has grown to ~ I5OM0 
after ~ 10^ yr, and to ~ TOOM© after ~ 3 x 10^ yr, the total lifetime of a very 
massive star. 



which is two orders of magnitude smaller than before. Several features arc 
evident in this plot. First, the central clump does not undergo further sub- 
fragmentation, and is likely to form a single Population III star. Second, a 
companion clump is visible at a distance of ^ 0.25 pc. If negative feedback 
from the first-forming star is ignored, this companion clump would undergo 
runaway collapse on its own approximately ~ 3 Myr later. This timescale is 
comparable to the lifetime of a very massive star (VMS) [59]. If the second 
clump was able to survive the intense radiative heating from its neighbor, it 
could become a star before the first one explodes as a supernova (SN). Whether 
more than one star can form in a low-mass halo thus crucially depends on the 
degree of synchronization of clump formation. Finally, the non-axisymmetric 
disturbance induced by the companion clump, as well as the angular momen- 
tum stored in the orbital motion of the binary system, allow the system to 
overcome the angular momentum barrier for the collapse of the central clump 
(see [216]). 

The recent discovery of stars like HE0107-5240 with a mass of 0.8Mq 
and an iron abundance of [Fe/H] = —5.3 [90] shows that at least some low 
mass stars could have formed out of extremely low-metallicity gas. The above 
simulations show that although the majority of clumps arc very massive, a 
few of them, like the secondary clump in Fig. 26, are significantly less massive. 
Alternatively, low-mass fragments could form in the dense, shock-compressed 
shells that surround the first hypernovae [234]. 
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How massive were the first stars? Star formation typically proceeds from 
the 'inside-out', through the accretion of gas onto a central hydrostatic core. 
Whereas the initial mass of the hydrostatic core is very similar for primordial 
and present-day star formation [274], the accretion process - ultimately re- 
sponsible for setting the final stellar mass, is expected to be rather different. 
On dimensional grounds, the accretion rate is simply related to the sound 
speed cubed over Newton's constant (or equivalently given by the ratio of the 
Jeans mass and the free-fall time): Mace ^ c^/G oc T'^/'^, A simple compari- 
son of the temperatures in present-day star forming regions (T ~ 10 K) with 
those in primordial ones (T ~ 200 — 300 K) already indicates a difference in 
the accretion rate of more than two orders of magnitude. 

The above refined simulation enables one to study the three-dimensional 
accretion flow around the protostar (see also [276, 304, 357]). The gas may 
now reach densities of 10^^ cm"'' before being incorporated into a central sink 
particle. At these high densities, three-body reactions [280] convert the gas 
into a fully molecular form. Figure 27 shows how the molecular core grows 
in mass over the first ^ 10"* yr after its formation. The accretion rate {left 
panel) is initially very high, M^cc ~ O.IM© yr~^, and subsequently declines 
according to a power law, with a possible break at ~ 5000 yr. The mass of the 
molecular core {right panel) ^ taken as an estimator of the proto-stellar mass, 
grows approximately as: ~ J Maccdi oc t^-^^ . A rough upper limit for the 
final mass of the star is then: M^{t = 3 x lO^yr) 7OOM0. In deriving this 
upper bound, we have conservatively assumed that accretion cannot go on for 
longer than the total lifetime of a massive star. 

Can a Population III star ever reach this asymptotic mass limit? The an- 
swer to this question is not yet known with any certainty, and it depends on 
whether the accretion from a dust-free envelope is eventually terminated by 
feedback from the star (e.g., [276, 304, 357, 277]). The standard mechanism 
by which accretion may be terminated in metal-rich gas, namely radiation 
pressure on dust grains [386] , is evidently not effective for gas with a primor- 
dial composition. Recently, it has been speculated that accretion could instead 
be turned off through the formation of an H II region [277], or through the 
photo-evaporation of the accretion disk [357]. The termination of the accre- 
tion process defines the current unsolved frontier in studies of Population III 
star formation. Current simulations indicate that the first stars were pre- 
dominantly very massive (> 30Mq), and consequently rather different from 
present-day stellar populations. The crucial question then arises: How and 
when did the transition take place from the early formation of massive stars 
to that of low-mass stars at later times? We address this problem next. 

The very first stars, marking the cosmic Renaissance of structure forma- 
tion, formed under conditions that were much simpler than the highly complex 
environment in present-day molecular clouds. Subsequently, however, the sit- 
uation rapidly became more complicated again due to the feedback from the 
first stars on the IGM. Supernova explosions dispersed the nucleosynthetic 
products from the first generation of stars into the surrounding gas (e.g., 
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[241, 261, 361]), including also dust grains produced in the explosion itself 
[222, 364]. Atomic and molecular cooling became much more efficient after 
the addition of these metals. Moreover, the presence of ionizing cosmic rays, 
as well as of UV and X-ray background photons, modified the thermal and 
chemical behavior of the gas in important ways (e.g., [232, 233]). 

Early metal enrichment was likely the dominant effect that brought about 
the transition from Population III to Population II star formation. Recent 
numerical simulations of collapsing primordial objects with overall masses of 
~ IO^Mq, have shown that the gas has to be enriched with heavy elements 
to a minimum level of Zait — W~^'^Zq, in order to have any effect on the 
dynamics and fragmentation properties of the system [275, 60, 64]. Normal, 
low-mass (Population II) stars are hypothesized to only form out of gas with 
metallicity Z > Zait- Thus, the characteristic mass scale for star formation 
is expected to be a fimction of metallicity, with a discontinuity at Zcrit where 
the mass scale changes by ~ two orders of magnitude. The redshift where this 
transition occurs has important implications for the early growth of cosmic 
structure, and the resulting observational signature (e.g., [392, 141, 234, 322]) 
include the extended nature of reionization [144]. 

For additional detailes about the properties of the first stars, see the com- 
prehensive review by Bromm & Larson (2004) [66]. 

5.3 Gamma-ray Bursts: Probing the First Stars One Star at a 
Time 

Gamma-Ray Bursts (GRBs) are believed to originate in compact remnants 
(neutron stars or black holes) of massive stars. Their high luminosities make 
them detectable out to the edge of the visible Universe [94, 214]. GRBs offer 
the opportunity to detect the most distant (and hence earliest) population of 
massive stars, the so-called Population III (or Pop III), one star at a time. In 
the hierarchical assembly process of halos which are dominated by cold dark 
matter (CDM), the first galaxies should have had lower masses (and lower 
stellar luminosities) than their low-redshift counterparts. Consequently, the 
characteristic luminosity of galaxies or quasars is expected to decline with 
increasing redshift. GRB afterglows, which already produce a peak flux com- 
parable to that of quasars or starburst galaxies at 2; ~ 1 — 2, are therefore 
expected to outshine any competing source at the highest redshifts, when the 
first dwarf galaxies have formed in the Universe. 

The first-year polarization data from the Wilkinson Microwave Anisotropy 
Probe ( WMAP) indicates an optical depth to electron scattering of ~ 17 ± 
4% after cosmological recombination [203, 348]. This implies that the first 
stars must have formed at a redshift z ~10— 20, and reionized a substantial 
fraction of the intergalactic hydrogen around that time [83, 93, 345, 394, 403]. 
Early reionization can be achieved with plausible star formation parameters 
in the standard ylCDM cosmology: in fac;t, the required optical depth can 
be achieved in a variety of very different ionization histories since WMAP 
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Fig. 28. Illustration of a long-duration gamma-ray burst in the popular "coUapsar" 
model. The collapse of the core of a massive star (which lost its hydrogen envelope) 
to a black hole generates two opposite jets moving out at a speed close to the speed 
of light. The jets drill a hole in the star and shine brightly towards an observer 
who happened to be located within with the coUimation cones of the jets. The jets 
emenating from a single massive star are so bright that they can be seen across the 
Universe out to the epoch when the first stars have formed. Upcoming observations 
by the Swift satellite will have the sensitivity to reveal whether the first stars served 
as progenitors of gamma-ray bursts (for updates see http://swift.gsfc.nasa.gov/). 

places only an integral constraint on these histories [176]. One would hke to 
probe the fuU history of reionization in order to disentangle the properties 
and formation history of the stars that are responsible for it. GRB afterglows 
offer the opportunity to detect stars as well as to probe the metal enrichment 
level [141] of the intervening IGM. 

GRBs, the electromagnetically-brightest explosions in the Universe, should 
be detectable out to redshifts z > 10 [94, 214]. High-redshift GRBs can be 
identified through infrared photometry, based on the Lya break induced by ab- 
sorption of their spectrum at wavelengths below 1.216 fim [(H-z)/10]. Follow- 
up spectroscopy of high-redshift candidates can then be performed on a IO- 
meter-class telescope. Recently, the ongoing Swift mission [147] has detected 
a GRB originating at z ~ 6.3 (e.g., [179]), thus demonstrating the viability of 
GRBs as probes of the early Universe. 

There are four main advantages of GRBs relative to traditional cosmic 
sources such as quasars: 

(i) The GRB afterglow flux at a given observed time lag after the 7-ray 
trigger is not expected to fade significantly with increasing redshift, since 
higher redshifts translate to earlier times in the source frame, during which the 
afterglow is intrinsically brighter [94] . For standard afterglow lightcurves and 
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spectra, the increase in the luminosity distance with redshift is compensated 
by this cosmological time- stretching effect. 
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Fig. 29. GRB afterglow flux as a function of time since the 7-ray trigger in the ob- 
server frame (taken from [67]). The flux (solid curves) is calculated at the redshifted 
Lya wavelength. The dotted curves show the planned detection threshold for the 
James Webb Space Telescope (JWST), assuming a spectral resolution R — 5000 with 
the near infrared spectrometer, a signal to noise ratio of 5 per spectral resolution 
element, and an exposure time equal to 20% of the time since the GRB explosion 
(see http://www.ngst.stsci.edu/nms/main/ ). Each set of curves shows a sequence 
of redshifts, namely z = 5, 7, 9, 11, 13, and 15, respectively, from top to bottom. 

(ii) As already mentioned, in the standard ylCDM cosmology, galaxies form 
hierarchically, starting from small masses and increasing their average mass 
with cosmic time. Hence, the characteristic mass of quasar black holes and the 
total stellar mass of a galaxy were smaller at higher redshifts, making these 
sources intrinsically fainter [391]. However, GRBs are believed to originate 
from a stellar mass progenitor and so the intrinsic luminosity of their engine 
should not depend on the mass of their host galaxy. GRB afterglows are 
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therefore expected to outshine their host galaxies by a factor that gets larger 
with increasing rcdshift. 

(Hi) Since the progenitors of GRBs are believed to be stellar, they likely orig- 
inate in the most common star-forming galaxies at a given redshift rather 
than in the most massive host galaxies, as is the case for bright quasars [26]. 
Low-mass host galaxies induce only a weak ionization effect on the surround- 
ing IGM and do not greatly perturb the Hubble flow around them. Hence, 
the Lya damping wing should be closer to the idealized unperturbed IGM 
case and its detailed spectral shape should be easier to interpret. Note also 
that unlike the case of a quasar, a GRB afterglow can itself ionize at most 
^ 4 X IO'^'E^iMq of hydrogen if its UV energy is Er,i in units of 10^^ ergs 
(based on the available number of ionizing photons), and so it should have a 
negligible cosmic eff'ect on the surrounding IGM. 

(iv) GRB afterglows have smooth (broken power-law) continuum spectra un- 
like quasars which show strong spectral features (such as broad emission lines 
or the so-called "blue bump" ) that complicate the extraction of IGM absorp- 
tion features. In particular, the contimmm extrapolation into the Lya damp- 
ing wing (the so-called Gunn-Peterson absorption trough) during the epoch 
of reionization is much more straightforward for the smooth UV spectra of 
GRB afterglows than for quasars with an underlying broad Lya emission line 
[26]. 

The optical depth of the uniform IGM to Lya absorption is given by (Gunn 
& Peterson 1965 [163]), 

Tre^/aAariHi (2;s) _ f nbh\ f nm\^^^^ f - ^^^"^ 

- — ~' 6.45 X 10 xni 



mecH{zs) ■ Vo-03y Vo.3y V 10 

(107) 

where H « lOO/i km s~^ Mpc~^i7m^(l + ^s)^/^ is the Hubble parameter at 
the source redshift z, » 1, /„ = 0.4162 and A^, = 1216A are the os- 
cillator strength and the wavelength of the Lya transition; nni {zs) is the 
neutral hydrogen density at the source redshift (assuming primordial abun- 
dances); ilm and f2b are the present-day density parameters of all matter 
and of baryons, respectively; and Xm is the average fraction of neutral hy- 
drogen. In the second equality we have implicitly considered high-redshifts, 
{1 + z) max [(1 — Qjn — ^A)/^m, (^yi/^^m)^^^] , at which the vacuum en- 
ergy density is negligible relative to matter (i7yi -C fim) and the Universe is 
nearly flat; for i?™ = 0.3, i7,4 = 0.7 this corresponds to the condition z^ 1.3 
which is well satisfied by the reionization redshift. 

At wavelengths longer than Lya at the source, the optical depth obtains 
a small value; these photons redshift away from the line center along its red 
wing and never resonate with the line core on their way to the observer. The 
red damping wing of the Gunn-Peterson trough (Miralda-Escude 1998 [254]) 
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Although the nature of the central engine that powers the relativistic jets 
of GRBs is still unknown, recent evidence indicates that long-duration GRBs 
trace the formation of massive stars (e.g., [365, 383, 45, 211, 47, 264]) and in 
particular that long-duration GRBs are associated with Type Ib/c supernovae 
[351]. Since the first stars in the Universe are predicted to be predominantly 
massive [5, 62, 66], their death might give rise to large numbers of GRBs at 
high redshifts. In contrast to quasars of comparable brightness, GRB after- 
glows are short-lived and release ~ 10 orders of magnitude less energy into the 
surrounding IGM. Beyond the scale of their host galaxy, they have a negligible 
effect on their cosmological environment^. Consequently, they are ideal probes 
of the IGM during the rcionization epoch. Their rest-frame UV spectra can be 
used to probe the ionization state of the IGM through the spectral shape of 
the Gunn-Peterson (Lya) absorption trough, or its metal enrichment history 
through the intersection of enriched bubbles of supernova (SN) ejecta from 
early galaxies [141]. Afterglows that are unusually bright (> lOmJy) at radio 
frequencies should also show a detectable forest of 21 cm absorption lines due 
to enhanced HI column densities in sheets, filaments, and collapsed minihalos 
within the IGM [76, 140]. 

Another advantage of GRB afterglows is that once they fade away, one 
may search for their host galaxies. Hence, GRBs may serve as signposts of 
the earliest dwarf galaxies that are otherwise too faint or rare on their own 
for a dedicated search to find them. Detection of metal absorption lines from 
the host galaxy in the afterglow spectrum, offers an unusual opportunity to 
study the physical conditions (temperature, metallicity, ionization state, and 
kinematics) in the interstellar medium of these high-redshift galaxies. As Fig- 
ure 30 indicates, damped Lya absorption within the host galaxy may mask 
the clear signature of the Gunn-Peterson trough in some galaxies [67]. A 
small fraction (~ 10) of the GRB afterglows are expected to originate at 
redshifts z > 5 [61, 68]. This subset of afterglows can be selected photomet- 
rically using a small telescope, based on the Lya break at a wavelength of 



^ Note, however, that feedback from a single GRB or supernova on the gas confined 
within early dwarf galaxies could be dramatic, since the binding energy of most 
galaxies at a > 10 is lower than lO^'^ ergs [23]. 
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Fig. 30. Expected spectral shape of the Lya absorption trough due to intergalactic 
absorption in GRB afterglows (taken from [67]). The spectrum is presented in terms 
of the flux density Fi, versus relative observed wavelength A\, for a source redshift 
z = 7 (assumed to be prior to the final reionization phase) and the typical halo mass 
M = 4 X 10* Mq expected for GRB host galaxies that cool via atomic transitions. 
Top panel: Two examples for the predicted spectrum including IGM HI absorption 
(both resonant and damping wing), for host galaxies with (i) an age ts ~ yr, 
a UV escape fraction /esc ~ 10% and a Scalo initial mass function (IMF) in solid 
curves, or (ii) ts ~ 10** yr, /esc = 90% and massive (> WOMq) Pop III stars in 
dashed curves. The observed time after the 7-ray trigger is one hour, one day, and 
ten days, from top to bottom, respectively. Bottom panel: Predicted spectra one 
day after a GRB for a host galaxy with ts = 10^ yr, /esc = 10% and a Scalo IMF. 
Shown is the unabsorbed GRB afterglow (dot-short dashed curve), the afterglow 
with resonant IGM absorption only (dot-long dashed curve), and the afterglow with 
full (resonant and damping wing) IGM absorption (solid curve). Also shown, with 
1.7 magnitudes of extinction, are the afterglow with full IGM absorption (dotted 
curve), and attempts to reproduce this profile with a damped Lya absorption system 
in the host galaxy (dashed curves). (Note, however, that damped absorption of this 
type could be suppressed by the ionizing effect of the afterglow UV radiation on 
the surrounding interstellar medium of its host galaxy [289].) Most importantly, the 
overall spectral shape of the Lya trough carries precious information about the 
neutral fraction of the IGM at the source redshift; averaging over an ensemble of 
sources with similar redshifts can reduce ambiguities in the interpretation of each 
case due to particular local effects. 
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1.216 yum [(1 + z)/10], caused by intergalactic HI absorption. The challenge in 
the upcoming years will be to follow-up on these candidates spectroscopically, 
using a large (10-meter class) telescope. GRB afterglows are likely to revolu- 
tionize observational cosmology and replace traditional sources like quasars, 
as probes of the IGM at z > 5. The near future promises to be exciting for 
GRB astronomy as well as for studies of the high-redshift Universe. 

It is of great importance to constrain the Pop III star formation mode, 
and in particular to determine down to which redshift it continues to be 
prominent. The extent of the Pop III star formation will affect models of the 
initial stages of reionization (e.g., [394, 93, 343, 403, 12]) and metal enrichment 
(e.g., [234, 141, 144, 320, 340]), and will determine whether planned surveys 
will be able to effectively probe Pop III stars (e.g., [319]). The constraints 
on Pop III star formation will also determine whether the first stars could 
have contributed a significant fraction to the cosmic near-IR backgroimd (e.g., 
[311, 310, 193, 242, 113]). To constrain high-redshift star formation from GRB 
observations, one has to address two major questions: 

(1) What is the signature of GRBs that originate in metal- free. Pop III pro- 
genitors? Simply knowing that a given GRB came from a high redshift is not 
sufficient to reach a definite conclusion as to the nature of the progenitor. 
Pregalactic metal enrichment was likely inhomogeneous, and we expect nor- 
mal Pop I and II stars to exist in galaxies that were already metal-enriched 
at these high redshifts [68]. Pop III and Pop I /II star formation is thus pre- 
dicted to have occurred concurrently at z > 5. How is the predicted high 
mass-scale for Pop HI stars reflected in the observational signature of the 
resulting GRBs? Preliminary results from numerical simulations of Pop HI 
star formation indicate that circumburst densities are systematically higher 
in Pop HI environments. GRB afterglows will then be much brighter than 
for conventional GRBs. In addition, due to the systematically increased pro- 
genitor masses, the Pop HI distribution may be biased toward long-duration 
events. 

(2) The modelling of Pop HI cosmic star formation histories has a number of 

free parameters, such as the star formation efficiency and the strength of the 
chemical feedback. The latter refers to the timescale for, and spatial extent 
of, the distribution of the first heavy elements that were produced inside of 
Pop HI stars, and subsequently dispersed into the IGM by supernova blast 
waves. Comparing with theoretical GRB redshift distributions one can use 
the GRB redshift distribution observed by Swift to calibrate the free model 
parameters. In particular, one can use this strategy to measure the redshift 
where Pop HI star formation terminates. 

Figures 31 and 32 illustrate these issues (based on [68]). Figure 32 leads 
to the robust expectation that ^ 10% of all Swift bursts should originate at 
z > 5. This prediction is based on the contribution from Population I/II stars 
which are known to exist even at these high redshifts. Additional GRBs could 
be triggered by Pop III stars, with a highly uncertain efficiency. Assuming 
that long-duration GRBs are produced by the coUapsar mechanism, a Pop III 
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Fig. 31. Theoretical prediction for the comoving star formation rate (SFR) in units 
of Mq yr~^ Mpc~^, as a function of redshift (from [68]). We assume that cooling 
in primordial gas is due to atomic hydrogen only, a star formation efficiency of 
?7* = 10%, and reionization beginning at Zrcion ~ 17. Solid line: Total comoving 
SFR. Dotted lines: Contribution to the total SFR from Pop I/II and Pop III for the 
case of weak chemical feedback. Dashed lines: Contribution to the total SFR from 
Pop I/II and Pop III for the case of strong chemical feedback. Pop III star formation 
is restricted to high redshifts, but extends over a significant range, Az ~ 10 — 15. 

star with a close binary companion provides a plausible GRB progenitor. The 
Pop III GRB efBciency, reflecting the probability of forming sufBciently close 
and massive binary systems, to lie between zero (if tight Pop III binaries do 
not exist) and ~ 10 times the empirically inferred value for Population I/II 
(due to the increased fraction of black hole forming progenitors among the 
massive Pop III stars). 

A key ingredient in determining the underlying star formation history 
from the observed GRB redshift distribution is the GRB luminosity function, 
which is only poorly constrained at present. The improved statistics provided 
by Swift will enable the construction of an empirical luminosity function. 
With an improved luminosity function it would be possible to re-calibrate the 
theoretical prediction in Figure 2 more reliably. 

In order to predict the observational signature of high-redshift GRBs, it is 
important to know the properties of the GRB host systems. Within variants 
of the popular CDM model for structure formation, where small objects form 
first and subsequently merge to build up more massive ones, the first stars are 
predicted to form at 2; 20-30 in minihalos of total mass (dark matter plus 
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Fig. 32. Predicted GRB rate to be observed by Swift (from [68]). Shown is the 
observed number of bursts per year, dNQ^s/dln{l + z), as a function of redshift. All 
rates are calculated with a constant GRB efficiency, ?7grb ~ 2 x 10~^ bursts Mq^, 
using the cosmic SFRs from Fig. 31. Dotted lines: Contribution to the observed 
GRB rate from Pop I/II aud Pop III for the case of weak chemical feedback. Dashed 
lines: Contribution to the GRB rate from Pop I/II and Pop III for the case of 
strong chemical feedback. Filled circle: GRB rate from Pop III stars if these were 
responsible for reionizing the Universe at « ~ 17. 



gas) ~ IO^Mq [359, 23, 403]. These objects are the sites for the formation of 
the first stars, and thus are the potential hosts of the highest-redshift GRBs. 
What is the environm,ent in which the earliest GRBs and their afterglows did 
occur? This problem breaks down into two related questions: (i) what type 
of stars (in terms of mass, metalUcity, and clustering properties) will form 
in each minihalo?, and (ii) how will the ionizing radiation from each star 
modify the density structure of the surrounding gas? These two questions are 
fundamentally intertwined. The ionizing photon production strongly depends 
on the stellar mass, which in turn is determined by how the accretion flow 
onto the growing protostar proceeds under the influence of this radiation field. 
In other words, the assembly of the Population III stars and the development 
of an HIT region around them proceed simultaneously, and affect each other. 
As a preliminary illustration, Figure 27 describes the photo-evaporation as a 
self-similar champagne flow [337] with parameters appropriate for the Pop III 
case. 

Notice that the central density is significantly reduced by the end of the 
life of a massive star, and that a central core has developed where the density 
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Fig. 33. Effect of photoheating from a Population III star on tiie density profile in 
a high-redshift minihalo (from [69]). The curves, labeled by the time after the onset 
of the central point source, are calculated according to a self-similar model for the 
expansion of an HIT region. Numerical simulations closely conform to this analytical 
behavior. Notice that the central density is significantly reduced by the end of the 
life of a massive star, and that a central core has developed where the density is 
constant. 

is iK^arly constant. Such a flat density profllc is markedly different from that 
created by stellar winds (p oc r~^). Winds, and consequently mass- loss, may 
not be important for massive Population III stars [18, 210], and such a flat 
density proflle may be characteristic of GRBs that originate from metal-free 
Population III progenitors. 

The flrst galaxies may be surrounded by a shell of highly enriched material 
that was carried out in a SN-drivcn wind (sec Fig. 34). A GRB in that galaxy 
may show strong absorption lines at a velocity separation associated with the 
wind velocity. Simulating these winds and calculating the absorption proflle in 
the featureless spectrum of a GRB afterglow, will allow us to use the observed 
spectra of high-z GRBs and directly probe the degree of metal enrichment in 
the vicinity of the first star forming regions (see [141] for a semi-analytic 
treatment). 

As the early afterglow radiation propagates through the interstellar en- 
vironment of the GRB, it will likely modify the gas properties close to the 
source; these changes could in turn be noticed as time-dependent spectral fea- 
tures in the spectrum of the afterglow and used to derive the properties of the 
gas cloud (density, metal abundance, and size). The UV afterglow radiation 
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Fig. 34. Supernova explosion in the high-redshift Universe (from [65]). The snapshot 
is taken ~ 10^ yr after the explosion with total energy .Esn — 10^^ ergs. We show 
the projected gas density within a box of linear size 1 kpc. The SN bubble has 
expanded to a radius of 200 pc, having evacuated most of the gas in the minihalo. 
Inset: Distribution of metals. The stellar ejecta {gray dots) trace the metals and are 
embedded in pristine metal-poor gas (black dots). 

can induce detectable changes to the interstellar absorption features of the 
host galaxy [289]; dust destruction could have occurred due to the GRB X- 
rays [375, 136], and molecules could have been destroyed near the GRB source 
[112]. Quantitatively, all of the effects mentioned above strongly depend on 
the exact properties of the gas in the host system. 

Most studies to date have assumed a constant efficiency of forming GRBs 
per unit mass of stars. This simplifying assumption could lead, under different 
circumstances, to an overestimation or an underestimation of the frequency of 
GRBs. Metal-free stars are thought to be massive [5, 62] and their extended 
envelopes may suppress the emergence of relativistic jets out of their surface 
(even if such jets are produced through the collapse of their core to a spinning 
black hole). On the other hand, low-metallicity stars are expected to have 
weak winds with little angular momentum loss during their evolution, and so 
they may preferentially yield rotating central configurations that make GRB 
jets after core collapse. 

What kind of metal-free, Pop III progenitor stars may lead to Long- 
duration GRBs appear to be associated with Type Ib/c supernovae [351], 
namely progenitor massive stars that have lost their hydrogen envelope. This 
requirement is explained theoretically in the coUapsar model, in which the 
relativistic jets produced by core collapse to a black hole are unable to emerge 
relativistically out of the stellar surface if the hydrogen envelope is retained 
[231]. The question then arises as to whether the lack of metal line-opacity 
that is essential for radiation-driven winds in metal-rich stars, would make a 
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Pop III star retain its hydrogen envelope, thus quenching any relativistic jets 
and GRBs. 

Aside from mass transfer in a binary system, individual Pop III stars could 
lose their hydrogen envelope due to either: (i) violent pulsations, particularly 
in the mass range 1OO-14OM0, or (ii) a wind driven by helium lines. The 
outer stellar layers are in a state where gravity only marginally exceeds radia- 
tion pressure due to electron-scattering (Thomson) opacity. Adding the small, 
but still non-ncgligiblc contribution from the bound-free opacity provided by 
singly- ionized helium, may be able to unbind the atmospheric gas. Therefore, 
mass-loss might occur even in the absence of dust or any heavy elements. 

5.4 Emission Spectrum of Metal- Free Stars 

The evolution of metal-free (Population III) stars is qualitatively different 
from that of enriched (Population I and II) stars. In the absence of the cata- 
lysts necessary for the operation of the CNO cycle, nuclear burning does not 
proceed in the standard way. At first, hydrogen burning can only occur via 
the inefiicient PP chain. To provide the necessary luminosity, the star has to 
reach very high central temperatures (T^ ~ 10*-^ K). These temperatures are 
high enough for the spontaneous turn-on of helium burning via the triplc-a 
process. After a brief initial period of triple-a burning, a trace amount of 
heavy elements forms. Subsequently, the star follows the CNO cycle. In con- 
structing main-sequence models, it is customary to assiimc that a trace mass 
fraction of metals {Z ~ 10~^) is already present in the star (El Eid 1983 [115]; 
Castellani et al. 1983 [78]). 

Figures 35 and 36 show the luminosity L vs. effective temperatm-c T for 
zero-age main sequence stars in the mass ranges of 2-9OM0 (Fig. 35) and 
100-IOOOMq (Fig. 36). Note that above ~ lOOMg the effective temperature 
is roughly constant, Teff ~ lO^K, implying that the spectrum is independent of 
the mass distribution of the stars in this regime (Bromm, Kudritzky, & Loeb 
2001 [59]). As is evident from these figures (see also Tumlinson & ShuU 2000 
[370]), both the effective temperature and the ionizing power of metal- free 
(Pop III) stars are substantially larger than those of metal-rich (Pop I) stars. 
Metal-free stars with masses > 2OM0 emit between 10^^ and 10^^ H I and He 
I ionizing photons per second per solar mass of stars, where the lower value 
applies to stars of ~ 2OM0 and the upper value applies to stars of > IOOM0 
(see Tumhnson & Shull 2000 [370] and Bromm et al. 2001 [59] for more details). 
Over a lifetime of ~ 3 x 10^ years these massive stars produce lO'^-lO''' ionizing 
photons per stellar baryon. However, this powerful UV emission is suppressed 
as soon as the interstellar medium out of which new stars form is enriched 
by trace amounts of metals. Even though the collapsed fraction of baryons is 
small at the epoch of reionization, it is likely that most of the stars responsible 
for the reionization of the Universe formed out of enriched gas. 
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Fig. 35. Luminosity vs. effective temperature for zero-age main sequences stars in 
the mass range of 2-9OM0 (from Tumlinson & Shull 2000 [370]). The curves show 
Pop I {Zq = 0.02) and Pop III stars of mass 2, 5, 8, 10, 15, 20, 25, 30, 35, 40, 50, 
60, 70, 80, and 90 Mq. The diamonds mark decades in mctallicity in the approach 
to Z = from 10"^ down to 10"^ at 2 Mq, down to 10"^" at 15 Mq, and down 
to 10"" at 90 Mq. The dashed line along the Pop III ZAMS assumes pure H-He 
composition, while the solid line (on the left) marks the upper MS with Zc = 10~^° 
for the M > 15 Mq models. Squares mark the points corresponding to pre-enriched 
evolutionary models from El Eid et al. (1983) [115] at 80 Mq and from Castellani 
et al. (1983) [78] for 25 Mq. 

Will it he possible to infer the initial mass function (IMF) of the first stars 
from spectroscopic observations of the first galaxies? Figure 37 compares the 
observed spectrum from a Salpeter IMF {dNi,/dM oc M~'^-^^) and a heavy 
IMF (with all stars more massive than IOOMq) for a galaxy at Zg = 10. 
The latter case follows from the assumption that each of the dense clumps 
in the simulations described in the previous section ends up as a single star 
with no significant fragmentation or mass loss. The difference between the 
plotted spectra cannot be confused with simple reddening due to normal dust. 
Another distinguishing feature of the IMF is the expected flux in the hydrogen 
and helium recombination lines, such as Lya and He II 1640 A, from the 
interstellar medium surrounding these stars. We discuss this next. 



5.5 Emission of Recombination Lines from the First Galcixies 

The hard UV emission from a star cluster or a quasar at high redshift is likely 
reprocessed by the surrounding interstellar medium, producing very strong 
recombination lines of hydrogen and helium (Oh 1999 [270]; Tumlinson & 
ShuU 2000 [370]; see also Baltz, Gnedin & Silk 1998 [17]). We define iVio„ 
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Fig. 36. Same as Figure 35 but for very massive stars above lOOM© (from Bromm, 
Kudritzki, & Loeb 2001 [59]). Left solid line: Pop III zero-age main sequence 
(ZAMS). Right solid line: Pop I ZAMS. In each case, stellar luminosity (in L©) 
is plotted vs. effective temperature (in K). Diamond-shaped symbols: Stellar masses 
along the sequence, from lOOM© (bottom) to IOOOM0 (top) in increments of lOOM© . 
The Pop III ZAMS is systematically shifted to higher effective temperature, with a 
value of ~ 10® K which is approximately independent of mass. The luminosities, on 
the other hand, are almost identical in the two cases. 

to be the production rate of ionizing photons by the source. The emitted 
luminosity L^™^ per unit stellar mass in a particular recombination line is 
then estimated to be 



where pj^™^ is the probability that a recombination leads to the emission of a 
photon in the corresponding line, v is the frequency of the line and Pcont ^^'^ 
Pune escape probabilities for the ionizing photons and the line photons, 

respectively. It is natural to assume that the stellar cluster is surrounded by a 
finite H II region, and hence that Pcont is close to zero [387, 302] . In addition, 
Piine is likely close to unity in the H II region, due to the lack of dust in the 
ambient metal- free gas. Although the emitted line photons may be scattered 
by neutral gas, they diffuse out to the observer and in the end survive if the 
gas is dust free. Thus, for simplicity, we adopt a value of unity for Py^^^- 

As a particular example we consider case B recombination which yields 
of about 0.65 and 0.47 for the Lya and He II 1640 A lines, respectively. 
These numbers correspond to an electron temperature of 3 x lO^K and 
an electron density of ~ 10^ — 10'^ cm~'^ inside the H II region [354]. For 
example, we consider the extreme and most favorable case of metal- free stars 
all of which are more massive than ^ IOOMq. In this case Lf"^^ = 1.7 x 10^'^ 
and 2.2 x 10'^^ erg s~^Mq ^ for the recombination luminosities of Lya and He 



-^line — Pline'i'^^^ionU - PcontjP] 



(111) 



First Light 71 




2 4 6 8 

Wavelength [firn] 

Fig. 37. Comparison of the predicted flux from a Pop III star cluster at Zs = 10 
for a Salpeter IMF (Tumlinson & Shull 2000 [370]) and a massive IMF (Bromm et 
al. 2001 [59]). Plotted is the observed flux (in njy per 10® M© of stars) vs. observed 
wavelength (in fim) for a flat Universe with Ha ~ 0.7 and h = 0.65. Solid line: The 
case of a heavy IMF. Dotted line: The flducial case of a standard Salpeter IMF. The 
cutoff below Xohs = 1216 A (1 + Zs) = 1.34^m is due to complete Gunn-Peterson 
absorption (which is artificially assumed to be sharp). Cleaxly, for the same total 
stellar mass, the observable flux is larger by an order of magnitude for stars which 
are biased towards having masses > lOOM© . 



II 1640 A per stellar mass [59]. A cluster of lO^M© in such stars would then 
produce 4.4 and 0.6 xlO^L© in the Lya and He II 1640 Alines. Comparably- 
high luminosities would be produced in other recombination lines at longer 
wavelengths, such as He II 4686 A and Hq [270, 271]. 

The rest-frame equivalent width of the above emission lines measured 
against the stellar continuum of the embedded star cluster at the line wave- 
lengths is given by 

(rem \ 
-f^j , (112) 

where L\ is the spectral luminosity per unit wavelength of the stars at the 
line resonance. The extreme case of metal-free stars which are more massive 
than IOOM0 yields a spectral luminosity per unit frequency L^, = 2.7 x 10^^ 
and 1.8 x lO^'^ erg }iz~^AlQ^ at the corresponding wavelengths [59]. Con- 
verting to Lx, this yields rest-frame equivalent widths of Wx = 3100 A and 
1100 A for Lya and He II 1640 A , respectively. These extreme emission equiv- 
alent widths are more than an order of magnitude larger than the expectation 
for a normal cluster of hot metal-free stars with the same total mass and a 
Salpeter IMF under the same assumptions concerning the escape probabili- 
ties and recombination [209]. The equivalent widths arc, of course, larger by a 
factor of (1 +Zs) in the observer frame. Extremely strong recombination lines. 
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such as Lya and He II 1640 A, are therefore expected to be an additional 

spectral signature that is unique to very massive stars in the early Universe. 
The strong recombination hues from the first luminous objects are potentially 
detectable with JWST[271]. 

6 Supermassive Black holes 
6.1 The Principle of Self- Regulation 

The fossil record in the present-day Universe indicates that every bulged 
galaxy hosts a supermassive black hole (BH) at its center [206]. This con- 
clusion is derived from a variety of techniques which probe the dynamics of 
stars and gas in galactic nuclei. The inferred BHs are dormant or faint most 
of the time, but ocassionally flash in a short burst of radiation that lasts for 
a small fraction of the Hubble time. The short duty cycle acounts for the fact 
that bright quasars arc much less abundant than their host galaxies, but it 
begs the more fundamental question: why is the quasar activity so brief? A 
natural explanation is that quasars are suicidal, namely the energy output 
from the BHs regulates their own growth. 

Supermassive BHs make up a small fraction, < 10"'^, of the total mass 
in their host galaxies, and so their direct dynamical impact is limited to the 
central star distribution where their gravitational influence dominates. Dy- 
namical friction on the background stars keeps the BH close to the center. 
Random fluctuations in the distribution of stars induces a Brownian motion 
of the BH. This motion can be decribed by the same Langcvin equation that 
captures the motion of a massive dust particle as it responds to random kicks 
from the much lighter molecules of air around it [86]. The characteristic speed 
by which the BH wanders around the center is small, ^ (m*/A'/BH)^^^o-*, 
where and Mbh are the masses of a single star and the BH, respectively, 
and (7* is the stellar velocity dispersion. Since the random force fluctuates 
on a dynamical time, the BH wanders across a region that is smaller by a 
factor of ~ (TO,t/MBH)^^^ than the region traversed by the stars inducing the 
fluctuating force on it. 

The dynamical insignificance of the BH on the global galactic scale is mis- 
leading. The gravitational binding energy per rest-mass energy of galaxies is 
of order ~ (u*/c)^ < 10~^. Since BH are relativistic objects, the gravitational 
binding energy of material that feeds them amounts to a substantial fraction 
its rest mass energy. Even if the BH mass occupies a fraction as small as 
~ 10~^ of the baryonic mass in a galaxy, and only a percent of the accreted 
rest-mass energy leaks into the gaseous environment of the BH, this slight 
leakage can unbind the entire gas reservoir of the host galaxy! This order- 
of-magnitude estimate explains why quasars are short lived. As soon as the 
central BH accretes large quantities of gas so as to significantly increase its 
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mass, it releases large amounts of energy that would suppress further accretion 
onto it. In short, the BH growth is self-regulated. 

The principle of self-regulation naturally leads to a correlation between 
the final BH mass, Mbh, and the depth of the gravitational potential well 
to which the surrounding gas is confined which can be characterized by the 
velocity dispersion of the associated stars, ~ a^. Indeed such a correlation is 
observed in the present-day Universe [368]. The observed power-law relation 
between Mbh and cr* can be generalized to a correlation between the BH mass 
and the circular velocity of the host halo, Vc [130], which in turn can be related 
to the halo mass, Mhaio, and redshift, z [394] 



where eo ~ 10 is a constant, and as before ( = [(i7m/]7^)(Z\c/187r^)], 
= [l-|-(/2^/r?„)(l + 0)-3]-i, Ac = 18'K^ + 82d-S9(f, and d= 12^-1. If 
quasars shine near their Eddington limit as suggested by observations of low 
and high-redshift quasars [134, 384], then the above value of implies that 
a fraction of ~ 5-10% of the energy released by the quasar over a galactic 
dynamical time needs to be captured in the surrounding galactic gas in order 
for the BH growth to be self-regulated [394] . 

With this interpretation, the Mbh-c* relation reflects the limit introduced 
to the BH mass by self-regulation; deviations from this relation are inevitable 
during episodes of BH growth or as a result of mergers of galaxies that have 
no cold gas in them. A physical scatter around this upper envelope could 
also result from variations in the efficiency by which the released BH energy 
couples to the surrounding gas. 

Various prescriptions for self-regulation were sketched by Silk & Rees [339] . 
These involve either energy or momentum-driven winds, where the latter type 
is a factor of ~ Vc/c less efficient [35, 199, 262]. Wyithe & Loeb [394] demon- 
strated that a particularly simple prescription for an energy-driven wind can 
reproduce the luminosity function of quasars out to highest measured red- 
shift, z ~ 6 (see Figs. 38 and 40), as well as the observed clustering properties 
of quasars at .2 ~ 3 [398] (see Fig. 41). The prescription postulates that: (i) 
self-regulation leads to the growth of Mbh up the redshift-independent limit 
as a function of Vc in Eq. (113), for all galaxies throughout their evolution; 
and (ii) the growth of A'/bh to the limiting mass in Eq. (113) occurs through 
halo merger episodes during which the BH shines at its Eddington luminosity 
(with the median quasar spectrum) over the dynamical time of its host galaxy, 
idyn- This model has only one adjustable parameter, namely the fraction of 
the released quasar energy that couples to the surrounding gas in the host 
galaxy. This parameter can be fixed based on the Mbh-c* relation in the lo- 
cal Univc;rsc [130]. It is remarkable that the combination of the above simple 
prescription and the standard ylCDM cosmology for the evolution and merger 
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Fig. 38. Comparison of the observed and model luminosity functions (from [394]). 
The data points at z < 4 are summarized in Ref. [286], while the light lines show 
the double power-law fit to the 2dF quasar luminosity function [56]. At z ^ 4.3 
and z ~ 6.0 the data is from Refs. [125]. The grey regions show the 1-a range 
of logarithmic slope ([—2.25, —3.75] at z ~ 4.3 and [—1.6, —3.1] at z ~ 6), and 
the vertical bars show the uncertainty in the normalization. The open circles show 
data points converted from the X-ray luminosity function [20] of low luminosity 
quasars using the median quasar spectral energy distribution. In each panel the 
vortical dashed lines correspond to the Eddington luminosities of BHs bracketing 
the observed range of the Mbh-Vc relation, and the vertical dotted line corresponds 
to a BH in a lO^-^M© galaxy. 

rate of galaxy halos, lead to a satisfactory agreement with the rich data set 
on quasar evolution over cosmic history. 

The cooling time of the heated gas is typically longer than its dynamical 
time and so the gas should expand into the galactic halo and escape the 
galaxy if its initial temperature exceeds the virial temperature of the galaxy 
[394]. The quasar remains active during the dynamical time of the initial 
gas reservoir, ^ 10^ years, and fades afterwards due to the dilution of this 
reservoir. Accretion is halted as soon as the quasar supplies the galactic gas 
with more than its binding energy. The BH growth may resume if the cold 
gas reservoir is replenished through a new merger. 

Following the early analytic work, extensive numerical simulations by 
Springel, Hernquist, & Di Matteo (2005) [350] (see also Di Matteo et al. 2005 
[108]) demonstrated that galaxy mergers do produce the observed correla- 
tions between black hole mass and spheroid properties when a similar energy 
feedback is incorporated. Because of the limited resolution near the galaxy 
nucleus, these simiilations adopt a simple prescription for the accretion flow 
that feeds the black hole. The actual feedback in reality may depend crucially 
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on the geometry of this flow and the physical mechanism that couples the 
energy or momentum output of the quasar to the surrounding gas. 



Fig. 39. Simulation images of a merger of galaxies resulting in quasar activity that 
eventually shuts-off the accretion of gas onto the black hole (from Di Matteo et 
al. 2005 [108]). The upper (lower) panels show a sequence of snapshots of the gas 
distribution during a merger with (without) feedback from a central black hole. The 
temperature of the gas is color coded. 

Agreement between the predicted and observed correlation function of 
quasars (Fig. 41) is obtained only if the BH mass scales with redshift as in 
Eq. (113) and the quasar lifetime is of the order of the dynamical time of the 
host galactic disk [398], 



This characterizes the timescale it takes low angular momentum gas to set- 
tle inwards and feed the black hole from across the galaxy before feedback 
sets in and suppresses additional infall. It also characterizes the timescale for 
establishing an outflow at the escape speed from the host spheroid. 

The inflow of cold gas towards galaxy centers during the growth phase 
of the BH would naturally be accompanied by a burst of star formation. 
The fraction of gas that is not consumed by stars or ejected by supernovae, 
will continue to feed the BH. It is therefore not surprising that quasar and 
starburst activities co-exist in Ultra Luminous Infrared Galaxies [356], and 
that all quasars show broad metal lines indicating a super-solar metallicity 
of the surrounding gas [106]. Applying a similar self- regulation principle to 
the stars, leads to the expectation [394, 197] that the ratio between the mass 
of the BH and the mass in stars is independent of halo mass (as observed 
locally [243]) but increases with redshift as oc ^(z)i/2(l + z)3/2. A consistent 
trend has indeed been inferred in an observed sample of gravitationally-lensed 
quasars [305]. 
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Fig. 40. The comoving density of supermassive BHs per unit BH mass (from [394]). 
The grey region shows the estimate based on the observed velocity distribution 
function of galaxies in Ref. [336] and the Mbh^c relation in Eq. (113). The lower 
bound corresponds to the lower limit in density for the observed velocity function 
while the grey lines show the extrapolation to lower densities. We also show the 
mass function computed at z = 1, 3 and 6 from the Press-Schecliter[292] halo mass 
function and Eq. (113), as well as the mass function at « ~ 2.35 and z ^ 3 implied 
by the observed density of quasars and a quasax lifetime of order the dynamical time 
of the host galactic disk, tdyn (dot-dashed lines). 

The upper mass of galaxies may also be regulated by the energy output 
from quasar activity. This would account for the fact that cooling flows are 
suppressed in present-day X-ray clusters [123, 91, 273], and that massive BHs 
and stars in galactic bulges were already formed at ^ 2. The quasars discov- 
ered by the Sloan Digital Sky Survey {SDSS) at 2 ~ 6 mark the early growth of 
the most massive BHs and galactic spheroids. The present-day abundance of 
galaxies capable of hosting BHs of mass ^ 10^ Mq (based on Eq. 113) already 
existed at 2 ~ 6 [225]. At some epoch, the quasar energy output may have 
led to the extinction of cold gas in these galaxies and the suppression of fur- 
ther star formation in them, leading to an apparent "anti-hierarchical" mode 
of galaxy formation where massive spheroids formed early and did not make 
new stars at late times. In the course of subsequent merger events, the cores of 
the most massive spheroids acquired an envelope of collisionless matter in the 
form of already-formed stars or dark matter [225], without the proportional 
accretion of cold gas into the central BH. The upper limit on the mass of the 
central BH and the mass of the spheroid is caused by the lack of cold gas and 
cooling flows in their X-ray halos. In the cores of cooling X-ray clusters, there 
is often an active central BH that supplies sufficient energy to compensate for 
the cooling of the gas [91, 123, 35]. The primary physical process by which 
this energy couples to the gas is still unknown. 
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Fig. 41. Predicted correlation function of quasars at various redshifts in comparison 
to the 2dF data [101] (from [398]). The dark lines show the correlation function 
predictions for quasars of various apparent B-band magnitudes. The 2dF limit is 
B ~ 20.85. The lower right panel shows data from entire 2dF sample in comparison 
to the theoretical prediction at the mean quasar redshift of (z) = 1.5. The B — 20.85 
prediction at this redshift is also shown by thick gray lines in the other panels to 
guide the eye. The predictions are based on the scaling Mbh oc in Eq. (113). 

6.2 Feedback on Large Intergalactic Scales 

Aside from affecting their host galaxy, quasars disturb their large-scale cos- 
mological environment. Powerful quasar outflows are observed in the form of 
radio jets [34] or broad-absorption-line winds [160]. The amount of energy 
carried by these outflows is largely unknown, but could be comparable to the 
radiative output from the same quasars. Furlanetto & Loeb [139] have calcu- 
lated the intergalactic volume filled by such outflows as a function of cosmic 
time (sec Fig. 42). This volume is likely to contain magnetic fields and metals, 
providing a natural source for the observed magnetization of the metal-rich 
gas in X-ray clusters [207] and in galaxies [103]. The injection of energy by 
quasar outflows may also explain the deficit of Lya absorption in the vicin- 
ity of Lyman-break galaxies [7, 100] and the required pre-heating in X-ray 
clusters [54, 91]. 

Beyond the reach of their outflows, the brightest SDSS quasars at ^ > 6 
are inferred to have ionized exceedingly large regions of gas (tens of comov- 
ing Mpc) around them prior to global reionization (see Fig. 43 and Refs. 
[381, 400]). Thus, quasars must have suppressed the faint-end of the galaxy 
luminosity function in these regions before the same occurred throughout the 
Universe. The recombination time is comparable to the Hubble time for the 
mean gas density at 2; ~ 7 and so ionized regions persist [272] on these large 
scales where inhomogeneities are small. The minimum galaxy mass is increased 
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Fig. 42. The global influence of magnetized quasar outflows on the intergalactic 
medium (from [139]). Upper Panel: Predicted volume filling fraction of magnetized 
quasar bubbles F{z), as a function of redshift. Lower Panel: Ratio of normalized 
magnetic energy density, UB/e-\, to the fiducial thermal energy density of the inter- 
galactic medium ufid — 3n{z)kTiGM , whore Tigm = 10'' K, as a function of redshift 
(see [139] for more details). In each panel, the solid curves assume that the blast 
wave created by quasar ouflows is nearly (80%) adiabatic, and that the minimum 
halo mass of galaxies, Mh,min, is determined by atomic cooling before reionization 
and by suppression due to galactic infall afterwards (top curve), Mh,min = 10^ Mq 
(middle curve), and Mh,min ~ W^'^Mq (bottom curve). The dashed curve assumes 
a fully-radiative blast wave and fixes Mh,min by the thresholds for atomic cooling 
and infall suppression. The vertical dotted line indicates the assumed redshift of 
complete reionization, Zr = 7. 



by at least an order of magnitude to a virial temperature of ~ lO^K in these 
ionized regions [23]. It would be particularly interesting to examine whether 
the faint end (cr^, < 30km s"'^) of the luminosity function of dwarf galaxies 
shows any moduluation on large-scales around rare massive BHs, such as M87. 

To find the volume filling fraction of relic regions from 2; ~ 6, we consider 
a BH of mass A^bh ^ 3 x 10^ Mq. Wc can estimate the comoving density of 
BHs directly from the observed quasar luminosity function and our estimate 
of quasar lifetime. At 2; ~ 6, quasars powered by Mbh ^ 3 x IO^Mq BHs had 
a comoving density of ^ 0.5Gpc~'^[394]. However, the Hubble time exceeds 
idyn by a factor of ~ 2 x 10^ (reflecting the square root of the density contrast 
of cores of galaxies relative to the mean density of the Universe), so that the 
comoving density of the bubbles created by the z ^ 6 BHs is ^ lO^Gpc"'^ (see 
Fig. 40). The density implies that the volume filling fraction of relic 2; ~ 6 
regions is small, < 10%, and that the nearest BH that had Mbh ^ 3 x lO^M© 
at 2; 6 (and could have been detected as an SDSS quasar then) should be 



First Light 79 



1 /3 

at a distance dbh ~ (47r/3 x 10^) Gpc ~ 140Mpc which is almost an order- 
of-magnitude larger than the distance of M87, a galaxy known to possess a 
BH of this mass [135]. 
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Fig. 43. Quasars serve as probes of the end of reionization. The measured size 
of the HIT regions around SDSS quasars can be used [396, 251] to demonstrate 
that a significant fraction of the intergalactic hydrogen was neutral at 2: ~ 6.3 
or else the inferred size of the quasar HIT regions would have been much larger 
than observed (assuming typical quasar lifetimes [248]). Also, quasars can be used 
to measure the redshift at which the intergalactic medium started to transmit Lya 
photons[381, 400]. The upper panel illustrates how the line-of-sight towards a quasar 
intersects this transition redshift. The resulting Lya transmission of the intrinsic 
quasar spectrum is shown schematically in the lower panel. 

What is the most massive BH that can be detected dynamically in a local 
galaxy redshift survey? SDSS probes a volume of ~ IGpc'^ out to a distance 
~ 30 times that of M87. At the peak of quasar activity at z ~ 3, the density of 
the brightest quasars implies that there should be ^ 100 BHs with masses of 
3x IO^^Mq per Gpc'^, the nearest of which will be at a distance dbh ~ 130Mpc, 
or ~ 7 times the distance to M87. The radius of gravitational influence of the 
BH scales as Mbh/i'c ^ ^'^hi^- We find that for the nearest 3 x lO^M© and 
3 X IO^'^Mq BHs, the angular radius of influence should be similar. Thus, the 
dynamical signature of ^ 3 x BHs on their stellar host should be 

detectable. 

6.3 What seeded the growth of the supermassive black holes? 

The BHs powering the bright quasars possess a mass of a few x lO^M©, 

and reside in galaxies with a velocity dispersion of ^ 500km s~^ [24]. A quasar 
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radiating at its Eddington limiting luminosity, Le — l-4x 10^^ erg s~^(A/bh/10®MQ), 
with a radiative efficiency, eiad = L^/Mc^ would grow exponentially in 
mass as a function of time i, Mbh = -^/gcod expji/i^} on a time scale, 
Ie = 4.1 X 10'' yr(erad/0.1). Thus, the required growth time in units of the 
Hubble time thubbio = 9 x 10^ yr[(l + z)/?]^^/^ is 

ihubbie 7 ) V^sood/lOOAfo; ^ ' 

The age of the Universe at z ~ 6 provides just sufficient time to grow an SDSS 
BH with Mbh ~ IO^Mq out of a stellar mass seed with erad = 10% [175]. The 
growth time is shorter for smaller radiative efficiencies, as expected if the seed 
originates from the optically-thick collapse of a supermassive star (in which 
case Msccd in the logarithmic factor is also larger). 
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Fig. 44. SPH simulation of the collapse of an early dwarf galaxy with a virial tem- 
perature just above the cooling threshold of atomic hydrogen and no H2 (from [63]). 
The image shows a snapshot of the gas density distribution at 2 ~ 10, indicating 
the formation of two compact objects near the center of the galaxy with masses of 
2.2 X lO^A^Q and 3.1 x lO'^M©, respectively, and radii < 1 pc. Sub-fragmentation 
into lower mass clumps is inhibited as long as molecular hydrogen is dissociated by 
a background UV flux. These circumstances lead to the formation of supermassive 
stars [220] that inevitably collapse and trigger the birth of supermassive black holes 
[220, 309]. The box size is 200 pc. 



What was the mass of the initial BH seeds? Were they planted in early 
dwarf galaxies through the collapse of massive, metal free (Pop-Ill) stars (lead- 
ing to Msccd of hundreds of solar masses) or through the collapse of even more 
massive, i.e. supermassive, stars [220] ? Bromm & Loeb [63] have shown 
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through a hydro dynamical simulation (see Fig. 44) that supermassive stars 
were likely to form in early galaxies at z 10 in which the virial temperature 
was close to the cooling threshold of atomic hydrogen, ~ lO^K. The gas in 
these galaxies condensed into massive ~ W^Mq clumps (the progenitors of 
supermassive stars), rather than fragmenting into many small clumps (the 
progenitors of stars), as it does in environments that are much hotter than 
the cooling threshold. This formation channel requires that a galax;y be close 
to its cooling threshold and immersed in a UV background that dissociates 
molecular hydrogen in it. These requirements should make this channel suffi- 
ciently rare, so as not to overproduce the cosmic mass density of supermassive 
BH. 

The minimum seed BH mass can be identified observationally through the 
detection of gravitational waves from BH binaries with Advanced LIGO [395] 
or with LISA [393]. Most of the mHz binary coalescence events originate at 
z > 1 ii the earliest galaxies included BHs that obey the Mbh-'Uc relation in 
Eq. (113). The number of LISA sources per unit redshift per year should drop 
substantially after reionization, when the minimum mass of galaxies increased 
due to photo-ionization heating of the intergalactic medium. Studies of the 
highest redshift sources among the few hundred detectable events per year, 
will provide unique information about the physics and history of BH growth 
in galaxies [327]. 

The early BH progenitors can also be detected as unresolved point sources, 
using the future James Webb Space Telescope {JWST). Unfortunately, the 
spectrum of metal-free massive and supermassive stars is the same, since their 
surface temperature ~ lO^K is independent of mass [59]. Hence, an unresolved 
cluster of massive early stars would show the same spectrum as a supermassive 
star of the same total mass. 

It is difficult to ignore the possible environmental impact of quasars on 
anthropic selection. One may wonder whether it is not a coincidence that our 
Milky- Way Galaxy has a relatively modest BH mass of only a few million solar 
masses in that the energy output from a much more massive (e.g. IO^Mq) 
black hole would have disrupted the evolution of life on our planet. A proper 
calculation remains to be done (as in the context of nearby Gamma- Ray Bursts 
[316]) in order to demonstrate any such link. 

7 Radiative Feedback from the First Sources of Light 
7.1 Escape of Ionizing Radiation from GalcLxies 

The intergalactic ionizing radiation field, a key ingredient in the development 
of reionization, is determined by the amount of ionizing radiation escaping 
from the host galaxies of stars and quasars. The value of the escape fraction 
as a function of redshift and galaxy mass remains a major uncertainty in all 
current studies, and could affect the cumulative radiation intensity by orders 
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of magnitude at any given redshift. Gas within halos is far denser than the 
typical density of the IGM, and in general each halo is itself embedded within 
an overdense region, so the transfer of the ionizing radiation must be followed 
in the densest regions in the Universe. Reionization simulations are limited in 
resolution and often treat the sources of ionizing radiation and their immediate 
surroundings as unresolved point sources within the large-scale intergalactic 
medium (see, e.g., Gnedin 2000 [152]). The escape fraction is highly sensitive 
to the three-dimensional distribution of the UV sources relative to the ge- 
ometry of the absorbing gas within the host galaxy (which may allow escape 
routes for photons along particular directions but not others). 

The escape of ionizing radiation {hiy > 13.6eV, A < 912 A) from the disks 
of present-day galaxies has been studied in recent years in the context of ex- 
plaining the extensive diffuse ionized gas layers observed above the disk in 
the Milky Way [300] and other galaxies [295, 183]. Theoretical models predict 
that of order 3-14% of the ionizing luminosity from O and B stars escapes the 
Milky Way disk [111, 110]. A similar escape fraction of /esc = 6% was deter- 
mined by Bland-Hawthorn & Maloney (1999) [46] based on Ha measurements 
of the Magellanic Stream. From Hopkins Ultraviolet Telescope observations of 
four nearby starburst galaxies (Leitherer et al. 1995 [217]; Hurwitz, Jelinsky, 
& Dixon 1997 [185]), the escape fraction was estimated to be in the range 
3%< /esc < 57%. If similar escape fractions characterize high redshift galax- 
ies, then stars could have provided a major fraction of the background radia- 
tion that reionized the IGM [236, 238]. However, the escape fraction from 
high-redshift galaxies, which formed when the Universe was much denser 
{p oc {1 + z)^), may be significantly lower than that predicted by models 
ment to describe present-day galaxies. Current reionization calculations as- 
sume that galaxies are isotropic point sources of ionizing radiation and adopt 
escape fractions in the range 5% < /esc < 60% [152]. 

Clumping is known to have a significant effect on the penetration and 
escape of radiation from an inhomogeneous medium [49, 385, 269, 173, 42]. 
The inclusion of dumpiness introduces several unknown parameters into the 
calculation, such as the number and ovcrdcnsity of the clumps, and the spa- 
tial correlation between the clumps and the ionizing sources. An additional 
complication may arise from hydrodynamic feedback, whereby part of the gas 
mass is expelled from the disk by stellar winds and supcrnovae (§8). 

Wood & Loeb (2000) [387] used a three-dimensional radiation transfer 
code to calculate the steady-state escape fraction of ionizing photons from 
disk galaxies as a function of redshift and galaxy mass. The gaseous disks were 
assumed to be isothermal, with a sound speed ~ 10 km s~^, and radially 
exponential, with a scale-length based on the characteristic spin parameter 
and virial radius of their host halos. The corresponding temperature of ^ 10^ 
K is typical for a gas which is continuousely heated by photo-ionization from 
stars. The sources of radiation were taken to be either stars embedded in the 
disk, or a central quasar. For stellar sources, the predicted increase in the 
disk density with redshift resulted in a strong decline of the escape fraction 
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with increasing redshift. The situation is different for a central quasar. Due 

to its higher luminosity and central location, the quasar tends to produce an 
ionization channel in the surrounding disk through which much of its ionizing 
radiation escapes from the host. In a steady state, only recombinations in this 
ionization channel must be balanced by ionizations, while for stars there are 
many ionization channels produced by individual star- forming regions and the 
total recombination rate in these channels is very high. Escape fractions > 10% 
were achieved for stars at z ^ 10 only if ^ 90% of the gas was expelled from 
the disks or if dense clumps removed the gas from the vast majority (> 80%) 
of the disk volume (see Fig. 45). This analysis applies only to halos with virial 
temperatures > 10** K. Ricotti & ShuU (2000) [302] reached similar conclusions 
but for a quasi-spherical configuration of stars and gas. They demonstrated 
that the escape fraction is substantially higher in low-mass halos with a virial 
temperature < 10'* K. However, the formation of stars in such halos depends 
on their uncertain ability to cool via the efficient production of molecular 
hydrogen. 
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Fig. 45. Escape fractions of stellar ionizing photons from a gaseous disk embedded 
within a lO^^M© halo which have formed at z = 10 (from Wood & Loeb 2000 
[387]). The curves show three different cases of dumpiness within the disk. The 
volume filling factor refers to either the ionizing emissivity, the gas clumps, or both, 
depending on the case. The escape fraction is substantial (> 1%) only if the gas 
distribution is highly clumped. 

The main uncertainty in the above predictions involves the distribution 
of the gas inside the host galaxy, as the gas is exposed to the radiation re- 
leased by stars and the mechanical energy deposited by supernovae. Given 
the fundamental role played by the escape fraction, it is desirable to calibrate 
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its value observationally. Steidel, Pettini, & Adelberger [352] reported a de- 
tection of significant Lyman continuum flux in the composite spectrum of 29 
Lyman break galaxies (LBG) with redshifts in the range z = 3.40 ± 0.09. 
They co-added the spectra of these galaxies in order to be able to measure 
the low flux. Another difficulty in the measurement comes from the need to 
separate the Lyman-limit break caused by the interstellar medium from that 
already produced in the stellar atmospheres. After correcting for intergalac- 
tic absorption, Steidel et al. [352] inferred a ratio between the emergent flux 
density at 1500A and 900A (rest frame) of 4.6 ± 1.0. Taking into account 
the fact that the stellar spectrum should already have an intrinsic Lyman 
discontinuity of a factor of ^ 3-5, but that only ^ 15 20% of the 1500 A pho- 
tons escape from typical LBGs without being absorbed by dust (Pettini et 
al. 1998 [290]; Adelberger et al. 2000 [6]), the inferred 900 A escape fraction 
is /esc ^ 10-20%. Although the galaxies in this sample were drawn from the 
bluest quartile of the LBG spectral energy distributions, the measurement im- 
plies that this quartile may itself dominate the hydrogen-ionizing background 
relative to quasars at ^ ~ 3. 

7.2 Propagation of Ionization Fronts in the IGM 

The radiation output from the first stars ionizes hydrogen in a growing volume, 
eventually encompassing almost the entire IGM within a single H II bubble. 
In the early stages of this process, each galaxy produces a distinct H II region, 
and only when the overall H II filling factor becomes significant do neighboring 
bubbles begin to overlap in large numbers, ushering in the "overlap phase" of 
reionization. Thus, the first goal of a model of reionization is to describe the 
initial stage, when each source produces an isolated expanding H II region. 

We assume a spherical ionized volume separated from the surrounding 
neutral gas by a sharp ionization front. Indeed, in the case of a stellar ioniz- 
ing spectrum, most ionizing photons are just above the hydrogen ionization 
threshold of 13.6 eV, where the absorption cross-section is high and a very 
thin layer of neutral hydrogen is sufficient to absorb all the ionizing photons. 
On the other hand, an ionizing source such as a quasar produces significant 
numbers of higher energy photons and results in a thicker transition region. 

In the absence of recombinations, each hydrogen atom in the IGM would 
only have to be ionized once, and the ionized proper volume would simply 
be determined by 

nuVp = , (116) 

where Uh is the mean number density of hydrogen and is the total number 
of ionizing photons produced by the source. However, the increased density 
of the IGM at high redshift implies that recombinations cannot be neglected. 
Indeed, in the case of a steady ionizing source (and neglecting the cosmological 
expansion), a steady-state volume would be reached corresponding to the 
Stromgren sphere, with recombinations balancing ionizations: 
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o^BulV, = ^ , (117) 

where the recombination rate depends on the square of the density and on 
the case B recombination coefficient as = 2.6 x 10~^^ cm'^ s~^ for hydrogen 
at T = 10^ K. The exact evolution for an expanding H II region, including 
a non-steady ionizing source, recombinations, and cosmological expansion, is 
given by (Shapiro & Giroux 1987 [329]) 

UH - iHv)l =^-ocB {n%) V, . (118) 

In this equation, the mean density uh varies with time as \/a?{t). A criti- 
cal physical ingredient is the dependence of recombination on the square of 
the density. This means that if the IGM is not uniform, but instead the gas 
which is being ionized is mostly distributed in high-density clumps, then the 
recombination time is very short. This is often dealt with by introducing a 
volume-averaged clumping factor C (in general time-dependent), defined by^ 

C={n%)/n'\j. (119) 

If the ionized volume is large compared to the typical scale of clumping, 

so that many clumps are averaged over, then equation (118) can be solved 
by supplementing it with equation (119) and specifying C. Switching to the 
comoving volume V, the resulting equation is 

dV 1 dN^ C n 

where the present number density of hydrogen is 

n% = 1.88 X 10-^ ( ^) cm-3 . (121) 

This number density is lower than the total number density of baryons n° by 
a factor of ^ 0.76, corresponding to the primordial mass fraction of hydrogen. 
The solution for Vi^) (generalized from Shapiro & Giroux 1987 [329]) around 
a source which turns on at t = ti is 

where 

F{t',t) = -aBn°,£^^dt" . (123) 



The recombination rate depends on the number density of electrons, and in using 
equation (119) we are neglecting the small contribution caused by partially or 
fully ionized helium. 
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At high redshift (when {1 + z) » \f2jn ^ — 1|), the scale factor varies as 



\ 2/3 



a(i)~ (^v^ifot) , (124) 



and with the additional assumption of a constant C the function F simplifies 
as follows. Defining 

m = ait)-^'^ , (125) 

we derive 

F{t',t) = ^ I-^^*') - = -0-262 [fit') - m] , (126) 



where the last equality assumes C = 10 and our standard choice of cosmo- 
logical parameters: i?^ = 0.3, J7/i ~ 0.7, and i?^ = 0.045. Although this 
expression for F{t',t) is in general an accurate approximation at high red- 
shift, in the particular case of the ylCDM model (where + fl^ = 1) we get 
the exact result by replacing equation (125) with 



The size of the resulting H II region depends on the halo which produces 
it. Consider a halo of total mass M and baryon fraction Qb/Qra- To derive 
a rough estimate, we assume that baryons are incorporated into stars with 
an efficiency of /star = 10%, and that the escape fraction for the resulting 
ionizing radiation is also /esc = 10%. If the stellar IMF is similar to the one 
measured locally [315], then N-y « 4000 ionizing photons are produced per 
baryon in stars (for a metallicity equal to 1/20 of the solar value). We define 
a parameter which gives the overall number of ionizations per baryon, 

.^ion = /star /esc • (128) 

If we neglect recombinations then we obtain the maximum comoving radius 
of the region which the halo of mass M can ionize, 

/3 N^V'"" /3 N,,^ My/3 (N,,^ M ^^'^ 

r-max = = \ir- ^n- -F, =680 kpc ' 



v47r n% J \Att n% Qrn mpj V 40 lO^M© 

(129) 

for our standard set of parameters. The actual radius never reaches this size if 
the recombination time is shorter than the lifetime of the ionizing source. For 
an instantaneous starburst with the Scalo (1998) [315] IMF, the production 
rate of ionizing photons can be approximated as 



dN-y _a-lN^ \ 1 if t < ts, 
dt a ts 1 ( 1 otherwise 



(130) 
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where = 4000, a = 4.5, and the most massive stars fade away with the 
characteristic timescale = 3x 10^ yr. In figure 46 we show the time evolution 
of the volume ionized by such a source, with the volume shown in units of 
the maximum volume Vmax which corresponds to Tmax in equation (129). We 
consider a source turning on at z = 10 (solid curves) or z = 15 (dashed curves), 
with three cases for each: no recombinations, C = 1, and C = 10, in order 
from top to bottom (Note that the result is independent of redshift in the case 
of no recombinations). When recombinations are included, the volume rises 
and reaches close to Vmax before dropping after the source turns off. At large t 
recombinations stop due to the dropping density, and the volume approaches 
a constant value (although V <C Vmax at large t ii C = 10). 




t [year] 

Fig. 46. Expanding H II region around an isolated ionizing source. The comov- 
ing ionized volume V is expressed in units of the maximum possible volume, 
Vmax = 47rrmax/3 [with Tmax giveu in equation (129)], and the time is measured 
after an instantaneous starburst which produces ionizing photons according to equa- 
tion (130). We consider a source turning on at 2 = 10 (solid curves) or 2 = 15 (dashed 
curves), with three cases for each: no recombinations, C = 1, and C = 10, in order 
from top to bottom. The no-recombination curve is identical for the different source 
redshifts. 

We obtain a similar result for the size of the H II region around a galaxy if 
we consider a mini-quasar rather than stars. For the typical quasar spectrum 
(Elvis et al. 1994 [122]), roughly 11,000 ionizing photons are produced per 
baryon incorporated into the black hole, assuming a radiative efficiency of 
6%. The efficiency of incorporating the baryons in a galaxy into a central 
black hole is low (< 0.6% in the local Universe, e.g. Magorrian et al. 1998 
[243]), but the escape fraction for quasars is likely to be close to unity, i.e., 
an order of magnitude higher than for stars (see previous sub-section). Thus, 
for every baryon in galaxies, up to 65 ionizing photons may be produced 
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by a central black hole and ~ 40 by stars, although both of these numbers 
for A^ion arc highly uncertain. These numbers suggest that in cither case the 
typical size of H II regions before reionization may be < 1 Mpc or ~ 10 Mpc, 
depending on whether IO^Mq halos or lO^^M© halos dominate. 

The ionization front around a bright transient source like a quasar expands 
at early times at nearly the speed of light. This occurs when the HII region 
is sufficiently small so that the production rate of ionizing photons by the 
central soiircc exceeds their consumption rate by hydrogen atoms within this 
volume. It is straightforward to do the accounting for these rates (including 
recombinations) taking the light propagation delay into account. This was 
done by Wyithc & Loeb [396] [see also White et al. (2003) [381]] who derived 
the general equation for the relativistic expansion of the comoving radius 
[r = (1 + z)rp] of the quasar H II region in an IGM with a neutral filling 
fraction xm (fixed by other ionizing sources) as, 



+ (1 + z) cxmnu 



(131) 



where c is the speed of light, C is the clumping factor, as = 2.6 x 10~^^cm^s~^ 
is the case-B recombination coefficient at the characteristic tcmpcratme of 
lO^K, and A'^ is the rate of ionizing photons crossing a shell at the radius 
of the HII region at time t. Indeed, for iV^ ^ oo the propagation speed of 
the proper radius of the HII region = r/(l + z) approaches the speed of 
light in the above expression, {drp/dt) c. The actual size of the HII region 
along the line-of-sight to a quasar can be inferred from the extent of the 
spectral gap between the quasar's rest-frame Lya wavelength and the start 
of Lya absorption by the IGM in the observed spectrum. Existing data from 
the SDSS quasars [396, 251, 401] provide typical values of Vp ~ 5Mpc and 
indicate for plausible choices of the quasar lifetimes that xhi > 0.1 at 2; > 6. 
These ionized bubbles could be imaged directly by future 21cm maps of the 
regions around the highest-redshift quasars [367, 397, 390] . 

The profile of the Lya emission line of galaxies has also been suggested as 
a probe of the ionization state of the IGM [223, 314, 81, 177, 240, 227, 246]. 
If the IGM is neutral, then the damping wing of the Gunn-Peterson trough in 
equation (108) is modified since Lya absorption starts only from the near edge 
of the ionized region along the line-of-sight to the source [81, 240]. Rhoads 
& Malhotra [246] showed that the observed abundance of galaxies with Lya 
emission at z ^ 6.5 indicates that a substantial fraction (tens of percent) 
of the IGM must be ionized in order to allow transmission of the observed 
Lya photons. However, if these galaxies reside in groups, then galaxies with 
peculiar velocities away from the observer will preferentially Doppler-shift the 
emitted Lya photons to the red wing of the Lya resonance and reduce the 
depression of the line profile [227, 85]. Additional uncertainties in the intrinsic 
line profile based on the geometry and the stellar or gaseous contents of the 
source galaxy [227, 314], as well as the clustering of galaxies which ionize 
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their immediate environment in groups [400, 145], limits this method from 
reaching robust conclusions. Imaging of the expected halos of scattered Lya 
radiation around galaxies embedded in a neutral IGM [223, 307] provide a 
more definitive test of the neutrality of the IGM, but is more challenging 
observationally. 

7.3 Reionization of Hydrogen 

In this section we summarize recent progress, both analytic and numerical, 
made toward elucidating the basic physics of reionization and the way in which 
the characteristics of reionization depend on the nature of the ionizing sources 
and on other input parameters of cosmological models. 

The process of the reionization of hydrogen involves several distinct stages. 
The initial, "pre-overlap" stage (using the terminology of Gnedin [152]) con- 
sists of individual ionizing sources turning on and ionizing their surroundings. 
The first galaxies form in the most massive halos at high redshift, and these 
halos are biased and are preferentially located in the highest-density regions. 
Thus the ionizing photons which escape from the galaxy itself (see §7.1) must 
then make their way through the surrounding high-density regions, which are 
characterized by a high recombination rate. Once they emerge, the ioniza- 
tion fronts propagate more easily into the low-density voids, leaving behind 
pockets of neutral, high-density gas. During this period the IGM is a two- 
phase medium characterized by highly ionized regions separated from neutral 
regions by ionization fronts. Furthermore, the ionizing intensity is very inho- 
mogeneous even within the ionized regions, with the intensity determined by 
the distance from the nearest source and by the ionizing luminosity of this 
source. 

The central, relatively rapid "overlap" phase of reionization begins when 
neighboring H II regions begin to overlap. Whenever two ionized bubbles 
are joined, each point inside their common boundary becomes exposed to 
ionizing photons from both sources. Therefore, the ionizing intensity inside 
H II regions rises rapidly, allowing those regions to expand into high-density 
gas which had previously recombined fast enough to remain neutral when 
the ionizing intensity had been low. Since each bubble coalescence accelerates 
the process of reionization, the overlap phase has the character of a phase 
transition and is expected to occur rapidly, over less than a Hubble time at 
the overlap redshift. By the end of this stage most regions in the IGM are 
able to see several unobscured sources, and therefore the ionizing intensity 
is much higher than before overlap and it is also much more homogeneous. 
An additional ingredient in the rapid overlap phase results from the fact that 
hierarchical structure formation models predict a galaxy formation rate that 
rises rapidly with time at the relevant redshift range. This process leads to 
a state in which the low-density IGM has been highly ionized and ionizing 
radiation reaches everywhere except for gas located inside self-shielded, high- 
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density clouds. This marks the end of the overlap phase, and this important 
landmark is most often referred to as the 'moment of rcionization'. 

Some neutral gas does, however, remain in high-density structures which 
correspond to Lyman Limit systems and damped Lya systems seen in ab- 
sorption at lower rcdshifts. The high-density regions are gradually ionized as 
galaxy formation proceeds, and the mean ionizing intensity also grows with 
time. The ionizing intensity continues to grow and to become more uniform 
as an increasing number of ionizing sources is visible to every point in the 
IGM. This "post-overlap" phase continues indefinitely, since collapsed objects 
retain neutral gas even in the present Universe. The IGM does, however, reach 
another milestone at z ^ 1.6, the breakthrough redshift [239]. Below this red- 
shift, all ionizing sources are visible to each other, while above this redshift 
absorption by the Lya forest implies that only sources in a small redshift 
range are visible to a typical point in the IGM. 

Semi-analytic models of the pre-overlap stage focus on the evolution of 
the H II filling factor, i.e., the fraction of the volume of the Universe which is 
filled by H II regions. We distinguish between the naive filling factor Fa u and 
the actual filling factor or porosity Qu n- The naive filling factor equals the 
number density of bubbles times the average volume of each, and it may exceed 
unity since when bubbles begin to overlap the overlapping volume is counted 
multiple times. However, as explained below, in the case of reionization the 
linearity of the physics means that Fu u is a very good approximation to 
Qh II up to the end of the overlap phase of rcionization. 

The model of individual H II regions presented in the previous section can 
be used to understand the development of the total filling factor. Starting with 
equation (120), if we assume a common clumping factor C for all H II regions 
then we can sum each term of the equation over all bubbles in a given large 
volume of the Universe, and then divide by this volume. Then V is replaced 
by the filling factor and A',^ by the total number of ionizing photons produced 
up to some time t, per unit volume. The latter quantity equals the mean 
number of ionizing photons per baryon times the mean density of baryons n^. 
Following the arguments leading to equation (129), we find that if we include 
only stars then 

^ = iVioni^col , (132) 
Ub 

where the collapse fraction Fcoi is the fraction of all the baryons in the Uni- 
verse which are in galaxies, i.e., the fraction of gas which settles into halos 
and cools efficiently inside them. In writing equation (132) we are assuming 
instantaneous production of photons, i.e., that the timescale for the forma- 
tion and evolution of the massive stars in a galaxy is short compared to the 
Hubble time at the formation redshift of the galaxy. In a model based on 
equation (120), the near-equality between Fu n and Qh u results from the 
linearity of this equation. First, the total number of ionizations equals the 
total number of ionizing photons produced by stars, i.e., all ionizing photons 
contribute regardless of the spatial distribution of sources; and second, the to- 
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tal recombination rate is proportional to the total ionized volume, regardless 
of its topology. Thus, even if two or more bubbles overlap the model remains 
an accurate approximation for Qh ii (at least until Qu u becomes nearly equal 
to 1). Note, however, that there still are a number of important simplifications 
in the model, including the assumption of a homogeneous (though possibly 
time-dependent) clumping factor, and the neglect of feedback whereby the 
formation of one galaxy may suppress further galaxy formation in neighbor- 
ing regions. These complications are discussed in detail below and in §7.5 and 
§8. 

Under these assumptions we convert equation (120), which describes indi- 
vidual H II regions, to an equation which statistically describes the transition 
from a neutral Universe to a fully ionized one (compare to Madau et al. 1999 
[239] and Haiman & Loeb 1997 [171]): 

^=076^-""^"-^""' ('^'^ 

where we assumed a primordial mass fraction of hydrogen of 0.76. The solution 
(in analogy with equation (122)) is 

where F{t',t) is determined by equations (123)-( 127). 

A simple estimate of the collapse fraction at high redshift is the mass 
fraction (given by equation (91) in the Press-Schechter model) in halos above 
the cooling threshold, which is the minimum mass of halos in which gas can 
cool efficiently. Assuming that only atomic cooling is effective during the red- 
shift range of reionization, the minimum mass corresponds roughly to a halo 
of virial temperature Tvir = 10^ K, which can be converted to a mass using 
equation (86). With this prescription we derive (for A'^ion = 40) the reioniza- 
tion history shown in Fig. 47 for the case of a constant clumping factor C. 
The solid curves show Qh ii as a function of redshift for a clumping factor 
C = (no recombinations), C = 1, C = 10, and C = 30, in order from 
left to right. Note that if C ~ 1 then recombinations are unimportant, but if 
C > 10 then recombinations significantly delay the reionization redshift (for 
a fixed star- formation history). The dashed curve shows the collapse fraction 
i^coi in this model. For comparison, the vertical dotted line shows the z = 5.8 
observational lower limit (Fan et al. 2000 [124]) on the reionization redshift. 

Clearly, star-forming galaxies in CDM hierarchical models are capable of 
ionizing the Universe at 2; ~ 6-15 with reasonable parameter choices. This 
has been shown by a large number of theoretical, semi-analytic calculations 
[138, 330, 171, 373, 89, 92, 392, 83, 371] as well as numerical simulations 
[79, 148, 152, 2, 296, 95, 342, 204, 186]. Similarly, if a small fraction (< 1%) 
of the gas in each galaxy accretes onto a central black hole, then the resulting 
mini-quasars are also able to reionize the Universe, as has also been shown 
using semi-analytic models [138, 172, 373, 392]. 
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Fig. 47. Semi-analytic calculation of the reionization of the IGM (for A'ion = 40), 
showing the redshift evolution of the filling factor Qh ii- Solid curves show Qh ii 
for a clumping factor C = (no recombinations), C = 1, C = 10, and C = 30, 
in order from left to right. The dashed curve shows the collapse fraction Fed, and 
the vertical dotted line shows the z = 5.8 observational lower limit (Fan et al. 2000 
[124]) on the reionization redshift. 

Although many models yield a reionization redshift around 7-12, the exact 
value depends on a number of uncertain parameters affecting both the source 
term and the recombination term in equation (133). The source parameters 
include the formation efficiency of stars and quasars and the escape fraction 
of ionizing photons produced by these sources. The formation efficiency of low 
mass galaxies may also be reduced by feedback from galactic outflows. These 
parameters affecting the sources are discussed elsewhere in this review (see 
§7.1, and 8). Even when the clumping is inhomogeneous, the recombination 
term in equation (133) is generally valid if C is defined as in equation (119), 
where we take a global volume average of the square of the density inside 
ionized regions (since neutral regions do not contribute to the recombination 
rate). The resulting mean clumping factor depends on the density and clus- 
tering of sources, and on the distribution and topology of density fluctuations 
in the IGM. Furthermore, the source halos should tend to form in overdense 
regions, and the clumping factor is affected by this cross-correlation between 
the sources and the IGM density. 

Miralda-Escude, Haehnelt, & Rees (2000) [256] presented a simple model 
for the distribution of density fluctuations, and more generally they discussed 
the implications of inhomogeneous clumping during reionization. They noted 
that as ionized regions grow, they more easily extend into low-density regions, 
and they tend to leave behind high-density concentrations, with these neutral 
islands being ionized only at a later stage. They therefore argued that, since 
at high-redshift the collapse fraction is low, most of the high-density regions. 
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which would dominate the clumping factor if they were ionized, will in fact 
remain neutral and occupy only a tiny fraction of the total volume. Thus, the 
development of reionization through the end of the overlap phase should occur 
almost exclusively in the low-density IGM, and the effective clumping factor 
during this time should be ^ 1, making recombinations relatively unimpor- 
tant (see Fig. 47). Only in the post-reionization phase, Miralda-Escude et al. 
(2000) [256] argued, do the high density clouds and filaments become gradu- 
ally ionized as the mean ionizing intensity further increases. 

The complexity of the process of reionization is illustrated by the numerical 
simulation of Gnedin [152] of stellar reionization (in ^CDM with = 0.3). 
This simulation uses a formulation of radiative transfer which relies on several 
rough approximations; although it does not include the effect of shadowing 
behind optically-thick clumps, it does include for each point in the IGM the 
effects of an estimated local optical depth around that point, plus a local 
optical depth around each ionizing source. This simulation helps to understand 
the advantages of the various theoretical approaches, while pointing to the 
complications which are not included in the simple models. Figures 48 and 
49, taken from Figure 3 in [152], show the state of the simulated Universe 
just before and just after the overlap phase, respectively. They show a thin 
(15 comoving kpc) slice through the box, which is 4 Mpc on a side, 
achieves a spatial resolution of lh~^ kpc, and uses 128^ each of dark matter 
particles and baryonic particles (with each baryonic particle having a mass 
of 5 X 10^ Mq). The figures show the redshift evolution of the mean ionizing 
intensity J21 (upper right panel), and visually the logarithm of the neutral 
hydrogen fraction (upper left panel), the gas density (lower left panel), and the 
gas temperature (lower right panel) . Note the obvious features resulting from 
the periodic boundary conditions assumed in the simulation. Also note that 
the intensity J21 is defined as the intensity at the Lyman limit, expressed 
in units of 10^^^ erg cm^^ s^^ sr~"'^Hz~"'^. For a given source emission, the 
intensity inside H II regions depends on absorption and radiative transfer 
through the IGM (e.g., Haardt & Madau 1996 [166]; Abel & Haehnelt 1999 
[1]) 

Figure 48 shows the two-phase IGM at 2: = 7.7, with ionized bubbles em- 
anating from one main concentration of sources (located at the right edge of 
the image, vertically near the center; note the periodic boundary conditions). 
The bubbles are shown expanding into low density regions and beginning to 
overlap at the center of the image. The topology of ionized regions is clearly 
complex: While the ionized regions are analogous to islands in an ocean of 
neutral hydrogen, the islands themselves contain small lakes of dense neutral 
gas. One aspect which has not been included in theoretical models of clump- 
ing is clear from the figure. The sources themselves arc located in the highest 
density regions (these being the sites where the earliest galaxies form) and 
must therefore ionize the gas in their immediate vicinity before the radiation 
can escape into the low density IGM. For this reason, the effective clump- 
ing factor is of order 100 in the simulation and also, by the overlap redshift. 
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Fig. 48. Visualization at z = 7.7 of a numerical simulation of reionization, adopted 
from Figure 3c of [152]. The panels display the logarithm of the neutral hydrogen 
fraction (upper left), the gas density (lower left), and the gas temperature (lower 
right). Also shown is the redshift evolution of the logarithm of the mean ionizing 
intensity (upper right). Note the periodic boundary conditions. 




Fig. 49. Visualization at z = 6.7 of a numerical simulation of reionization, adopted 
from Figure 3e of [152]. The panels display the logarithm of the neutral hydrogen 
fraction (upper left), the gas density (lower left), and the gas temperature (lower 
right). Also shown is the redshift evolution of the logarithm of the mean ionizing 
intensity (upper right). Note the periodic boundary conditions. 
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roughly ten ionizing photons have been produced per baryon. Figure 49 shows 

that by z = 6.7 the low density regions have all become highly ionized along 
with a rapid increase in the ionizing intensity. The only neutral islands left 
are the highest density regions (compare the two panels on the left). However, 
wc emphasize that the quantitative results of this simulation must be consid- 
ered preliminary, since the effects of increased resolution and a more accurate 
treatment of radiative transfer are yet to be explored. Methods are being de- 
veloped for incorporating a more complete treatment of radiative transfer into 
three dimensional cosmological simulations (e.g., [2, 296, 95, 342, 204, 186]). 

Gnedin, Ferrara, & Zweibel (2000) [151] investigated an additional effect 
of reionization. They showed that the Biermann battery in cosmological ion- 
ization fronts inevitably generates coherent magnetic fields of an amplitude 
~ 10~^^ Gauss. These fields form as a result of the breakout of the ionization 
fronts from galaxies and their propagation through the H I filaments in the 
IGM. Although the fields are too small to directly affect galaxy formation, 
they could be the seeds for the magnetic fields observed in galaxies and X-ray 
clusters today. 

If quasars contribute substantially to the ionizing intensity during reion- 
ization then several aspects of reionization are modified compared to the case 
of pure stellar reionization. First, the ionizing radiation emanates from a sin- 
gle, bright point-source inside each host galaxy, and can establish an escape 
route (H II funnel) more easily than in the case of stars which are smoothly 
distributed throughout the galaxy (§7.1). Second, the hard photons produced 
by a quasar penetrate deeper into the surrounding neutral gas, yielding a 
thicker ionization front. Finally, the quasar X-rays catalyze the formation of 

molecules and allow stars to keep forming in very small halos. 

Oh (1999) [270] showed that star-forming regions may also produce signif- 
icant X-rays at high redshift. The emission is due to inverse Compton scat- 
tering of CMB photons off rclativistic electrons in the ejecta, as well as ther- 
mal emission by the hot supernova remnant. The spectrum expected from 
this process is even harder than for typical quasars, and the hard photons 
photoionize the IGM efficiently by repeated secondary ionizations. The radia- 
tion, characterized by roughly equal energy per logarithmic frequency interval, 
would produce a uniform ionizing intensity and lead to gradual ionization and 
heating of the entire IGM. Thus, if this source of emission is indeed effective 
at high redshift, it may have a crucial impact in changing the topology of 
reionization. Even if stars dominate the emission, the hardness of the ionizing 
spectrum depends on the initial mass function. At high redshift it may be 
biased toward massive, efficiently ionizing stars, but this remains very much 
uncertain. 

Semi-analytic as well as numerical models of reionization depend on an 
extrapolation of hierarchical models to higher redshifts and lower-mass halos 
than the regime where the models have been compared to observations (see 
e.g. [392, 83, 371]). These models have the advantage that they are based on 
the current GDM paradigm which is supported by a variety of observations 
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of large-scale structure, galaxy clustering, and the CMB. The disadvantage is 
that the properties of high-rcdshift galaxies are derived from those of their 
host halos by prescriptions which are based on low redshift observations, and 
these prescriptions will only be tested once abundant data is available on 
galaxies which formed during the reionization era (see [392] for the sensitivity 
of the results to model parameters) . An alternative approach to analyzing the 
possible ionizing sources which brought about reionization is to extrapolate 
from the observed populations of galaxies and quasars at currently accessible 
redshifts. This has been attempted, e.g., by Madau et al. (1999) [239] and 
Miralda-Escude et al. (2000) [256]. The general conclusion is that a high- 
redshift source population similar to the one observed at z = 3 -4 would pro- 
duce roughly the needed ionizing intensity for reionization. However, Dijkstra, 
Haiman, & Loeb (2004) [107] constrained the role of quasars in reionizing the 
Universe based on the unresolved flux of the X-ray background. At any event, 
a precise conclusion remains elusive because of the same kinds of uncertain- 
ties as those found in the models based on CDM: The typical escape fraction, 
and the faint end of the luminosity function, are both not well determined 
even a,t z = 3-4, and in addition the clumping factor at high redshift must be 
known in order to determine the importance of recombinations. Future direct 
observations of the source population at redshifts approaching reionization 
may help resolve some of these questions. 

7.4 Photo-evaporation of Gaseous Halos After Reionization 

The end of the reionization phase transition resulted in the emergence of an 
intense UV background that filled the Universe and heated the IGM to tem- 
peratures of ~ 1-2 X lO^'K (see the previous section). After ionizing the rarefied 
IGM in the voids and filaments on large scales, the cosmic UV background 
penetrated the denser regions associated with the virialized gaseous halos of 
the first generation of objects. A major fraction of the collapsed gas had been 
incorporated by that time into halos with a virial temperature < lO^K, where 
the lack of atomic cooling prevented the formation of galactic disks and stars 
or quasars. Photoionization heating by the cosmic UV background could then 
evaporate much of this gas back into the IGM. The photo-evaporating halos, 
as well as those halos which did retain their gas, may have had a number of 
important consequences just after reionization as well as at lower redshifts. 

In this section we focus on the process by which gas that had already 
settled into virialized halos by the time of reionization was evaporated back 
into the IGM due to the cosmic UV background. This process was inves- 
tigated by Barkana fc Loeb (1999) [22] using semi-analytic methods and 
idealized numerical calculations. They first considered an isolated spherical, 
centrally-concentrated dark matter halo containing gas. Since most of the 
photo-evaporation occurs at the end of overlap, when the ionizing intensity 
builds up almost instantaneously, a sudden ilhunination by an external ion- 
izing background may be assumed. Self-shielding of the gas implies that the 
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halo interior sees a reduced intensity and a harder spectrum, since the outer 

gas layers preferentially block photons with energies just above the Lyman 
limit. It is useful to parameterize the external radiation field by a specific 
intensity per unit frequency, u, 



where Uj^ is the Lyman limit frequency, and J21 is the intensity at ex- 
pressed in units of 10~^^ erg cm~^ s~^ sr~^Hz~^. The intensity is normalized 
to an expected post-reionization value of around unity for the ratio of ioniz- 
ing photon density to the baryon density. Different power laws can be used to 
represent either quasar spectra (a ^ 1.8) or stellar spectra {a ~ 5). 

Once the gas is heated throughout the halo, some fraction of it acquires 
a sufficiently high temperature that it becomes unbound. This gas expands 
due to the resulting pressure gradient and eventually evaporates back into the 
IGM. The pressure gradient force (per unit volume) kS/ {T p / fj,mp) competes 
with the gravitational force of p GM/r"^ . Due to the density gradient, the ratio 
between the pressure force and the gravitational force is roughly equal to the 
ratio between the thermal energy ~ kT and the gravitational binding energy 
~ jjLViipGM /r (which is ^ fcTvir at the virial radius Tvir) per particle. Thus, 
if the kinetic energy exceeds the potential energy (or roughly if T > T^ir), 
the repulsive pressure gradient force exceeds the attractive gravitational force 
and expels the gas on a dynamical time (or faster for halos with T ^ Tvir). 

The left panel of Figure 50 (adopted from Fig. 3 of Barkana & Loeb 1999 
[22]) shows the fraction of gas within the virial radius which becomes unbound 
after reionization, as a function of the total halo circular velocity, with halo 
masses at z = 8 indicated at the top. The two pairs of curves correspond to 
spectral index a = 5 (solid) or a = 1.8 (dashed). In each pair, a calculation 
which assumes an optically-thin halo leads to the upper curve, but including 
radiative transfer and self-shielding modifies the result to the one shown by 
the lower curve. In each case self-shielding lowers the unbound fraction, but 
it mostly affects only a neutral core containing ~ 30% of the gas. Since high 
energy photons above the Lyman limit penetrate deep into the halo and heat 
the gas efficiently, a flattening of the spectral slope from a = 5 to a = 1.8 raises 
the unbound gas fraction. This figure is essentially independent of redshift if 
plotted in terms of circular velocity, but the conversion to a corresponding 
mass does vary with redshift. The characteristic circular velocity where most 
of the gas is lost is ^ 10 -15 km s^^, but clearly the effect of photo-evaporation 
is gradual, going from total gas removal down to no effect over a range of a 
factor of ~ 100 in halo mass. 

Given the values of the unbound gas fraction in halos of different masses, 
the Press-Schechter mass function (§4.1) can be used to calculate the total 
fraction of the IGM which goes through the process of accreting onto a halo 
and then being recycled into the IGM at reionization. The low-mass cutoff 
in this sum over halos is given by the lowest mass halo in which gas has 
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Fig. 50. Effect of photo-evaporation on individual halos and on the overall halo 
population. The left panel shows the unbound gas fraction (within the virial radius) 
versus total halo velocity dispersion or mass, adopted from Figure 3 of Barkana 
& Loeb (1999) [22]. The two pairs of curves correspond to spectral index a = 5 
(solid) or a = 1.8 (dashed), in each case at z = 8. In each pair, assuming an 
optically-thin halo leads to the upper curve, while the lower curve shows the result 
of including radiative transfer and self shielding. The right panel shows the total 
fraction of gas in the Universe which evaporates from halos at reionization, versus 
the reionization redshift, adopted from Figure 7 of Barkana & Loeb (1999) [22]. The 
solid line assumes a spectral index a — 1.8, and the dotted line assumes a = 5. 



assembled by the reionization redshift. This mass can be estimated by the 
linear Jeans mass Afj in equation (62). The Jeans mass does not in general 
precisely equal the limiting mass for accretion (see the discussion in the next 
section). Indeed, at a given redshift some gas can continue to fall into halos of 
lower mass than the Jeans mass at that redshift. On the other hand, the larger 
Jeans mass at higher redshifts means that a time-averaged Jeans mass may 
be more appropriate, as indicated by the filtering mass. In practice, the Jeans 
mass is sufficiently accurate since at z ~ 10-20 it agrees well with the values 
found in the numerical spherical collapse calculations of Haiman, Thoul, & 
Loeb (1996) [168]. 

The right panel of Figure 50 (adopted from Fig. 7 of Barkana & Loeb 1999 
[22]) shows the total fraction of gas in the Universe which evaporates from 
halos at reionization, versus the reionization redshift. The solid line assumes 
a spectral index a — 1.8, and the dotted line assumes a = 5, showing that 
the result is insensitive to the spectrum. Even at high redshift, the amount 
of gas which participates in photo-evaporation is significant, which suggests 
a number of possible implications as discussed below. The gas fraction shown 
in the figure represents most {^-^ 60-80% depending on the redshift) of the 
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collapsed fraction before reionization, although some gas does remain in more 
massive hales. 

The photo-evaporation of gas out of large numbers of halos may have 
interesting implications. First, gas which falls into halos and is expelled at 
reionization attains a different entropy than if it had stayed in the low-density 
IGM. The resulting overall reduction in the entropy is expected to be small - 
the same as would be produced by reducing the temperature of the entire IGM 
by a factor of ^ 1.5 - but localized effects near photo-evaporating halos may 
be more significant. Furthermore, the resulting ~ 20 km s~^ outflows induce 
small-scale fluctuations in peculiar velocity and temperature. These outflows 
are usually well below the resolution limit of most numerical simulations, but 
some outflows were resolved in the simulation of Bryan et al. (1998) [70]. The 
evaporating halos may consume a signiflcant number of ionizing photons in 
the post-overlap stage of reionization [174, 186], but a definitive determination 
requires detailed simulations which include the three-dimensional geometry of 
source halos and sink halos. 

Although gas is quickly expelled out of the smallest halos, photo-evaporation 
occurs more gradually in larger halos which retain some of their gas. These 
surviving halos initially expand but they continue to accrete dark matter and 
to merge with other halos. These evaporating gas halos could contribute to 
the high column density end of the Lya forest [51]. Abel & Mo (1998) [3] 
suggested that, based on the expected number of surviving halos, a large frac- 
tion of the Lyman limit systems at z ^ 3 may correspond to mini-halos that 
survived reionization. Surviving halos may even have identifiable remnants in 
the present Universe. These ideas thus off'er the possibility that a population 
of halos which originally formed prior to reionization may correspond almost 
directly to several populations that are observed much later in the history 
of the Universe. However, the detailed dynamics of photo-evaporating halos 
are complex, and detailed simulations arc required to confirm these ideas. 
Photo-evaporation of a gas cloud has been followed in a two dimensional sim- 
ulation with radiative transfer, by Shapiro & Raga (2000) [331]. They found 
that an evaporating halo would indeed appear in absorption as a damped 
Lya system initially, and as a weaker absorption system subsequently. Future 
simulations [186] will clarify the contribution to quasar absorption lines of the 
entire population of photo-evaporating halos. 

7.5 Suppression of the Formation of Low Mass Galaxies 

At the end of overlap, the cosmic ionizing background increased sharply, and 
the IGM was heated by the ionizing radiation to a temperature > 10^ K. Due 
to the substantial increase in the IGM temperature, the intergalactic Jeans 
mass increased dramatically, changing the minimum mass of forming galaxies 
[299, 117, 148, 255]. 

Gas infall depends sensitively on the Joans mass. When a halo more mas- 
sive than the Jeans mass begins to form, the gravity of its dark matter over- 
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comes the gas pressure. Even in halos below the Jeans mass, although the 
gas is initially held up by pressure, once the dark matter collapses its in- 
creased gravity pulls in some gas [168]. Thus, the Jeans mass is generally 
higher than the actual limiting mass for accretion. Before reionization, the 
IGM is cold and neutral, and the Jeans mass plays a secondary role in limit- 
ing galaxy formation compared to cooling. After reionization, the Jeans mass 
is increased by several orders of magnitude due to the photoionization heat- 
ing of the IGM, and hence begins to play a dominant role in limiting the 
formation of stars. Gas infall in a reionized and heated Universe has been 
investigated in a number of numerical simulations. Thoul & Weinberg (1996) 
[363] inferred, based on a spherically-symmetric collapse simulation, a reduc- 
tion of ~ 50% in the collapsed gas mass due to heating, for a halo of circular 
velocity ~ 50 km s~^ at 2; = 2, and a complete suppression of infall below 
Vc ^ 30 km s""'^. Kitayama & Ikeuchi (2000) [201] also performed spherically- 
symmetric simulations but included self-shielding of the gas, and found that 
it lowers the circular velocity thresholds by ~ 5 km s~^. Three dimensional 
numerical simulations [294, 378, 267] found a significant suppression of gas 
infall in even larger halos (14 ~ 75 km s~^), but this was mostly due to a 
suppression of late infall at 2; < 2. 

When a volume of the IGM is ionized by stars, the gas is heated to a 
temperature Tjgm ~ 10** K. If quasars dominate the UV background at reion- 
ization, their harder photon spectrum leads to Tjgm > 2 x 10^ K. Including 
the effects of dark matter, a given temperature results in a linear Jeans mass 
corresponding to a halo circular velocity of 



where we used equation (85) and assumed /i = 0.6. In halos with Vc > Vj, 
the gas fraction in infalling gas equals the universal mean of f2h/f2m, but 
gas infall is suppressed in smaller halos. Even for a small dark matter halo, 
once it collapses to a virial overdensity of Ac/fi^ relative to the mean, it 
can pull in additional gas. A simple estimate of the limiting circular velocity, 
below which halos have essentially no gas infall, is obtained by substituting 
the virial overdensity for the mean density in the definition of the Jeans mass. 
The resulting estimate is 



This value is in rough agreement with the numerical simulations mentioned 
before. A more recent study by Dijkstra et al. (2004) [107] indicates that at 
the high rcdshifts of z > 10 gas could nevertheless assemble into halos with 
circular velocities as low as Vc ~ 10 km s~^, even in the presence of a UV 
background. 
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Although the Jeans mass is closely related to the rate of gas infall at a 
given time, it does not directly yield the total gas residing in halos at a given 
time. The latter quantity depends on the entire history of gas accretion onto 
halos, as well as on the merger histories of halos, and an accurate description 
must involve a time-avcragcd Jeans mass. Gnedin [153] showed that the gas 
content of halos in simulations is well fit by an expression which depends 
on the filtering mass, a particular time-averaged Jeans mass (Gnedin fc Hui 
1998 [150]). Gnedin [153] calculated the Jeans and filtering masses using the 
mean temperature in the simulation to define the sound speed, and found the 
following fit to the simulation results: 

' [l+(2i/3-l)Mc/M]" 

where Mg is the average gas mass of all objects with a total mass M, fh = 
Qb/Qm is the universal baryon fraction, and the characteristic mass Mc is 
the total mass of objects which on average retain 50% of their gas mass. The 
characteristic mass was well fit by the filtering mass at a range of redshifts 
from 2; = 4 up to 2; ~ 15. 

The reionization process was not perfectly synchronized throughout the 
Universe. Large-scale regions with a higher density than the mean tend to form 
galaxies first and reionize earlier than underdense regions (see detailed discus- 
sion in §168). The suppression of low-mass galaxies by reionization will there- 
fore be modulated by the fluctuations in the timing of reionization. Babich & 
Loeb (2005) [14] considered the effect of inhomogeneous reionization on the 
power-spectrum of low-mass galaxies. They showed that the shape of the high 
redshift galaxy power spectrum on small scales in a manner which depends 
on the details of epoch of reionization. This effect is significantly larger than 
changes in the galaxy power spectrum due to the current uncertainty in the 
inflationary parameters, such as the tilt of the scalar power spectrum n and 
the running of the tilt a. Therefore, future high redshift galaxies surveys hop- 
ing to constrain inflationary parameters must properly model the effects of 
reionization, but conversely they will also be sensitive to the thermal history 
of the high redshift intergalactic medium. 



8 Feedback from Galactic Outflows 

8.1 Propagation of Supernova Outflow^s in the IGM 

Star formation is accompanied by the violent death of massive stars in su- 
pernova explosions. In general, if each halo has a fixed baryon fraction and a 
fixed fraction of the baryons turns into massive stars, then the total energy in 
supernovae outflows is proportional to the halo mass. The binding energy of 
the gas in the halo is proportional to the halo mass squared. Thus, outflows 
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are expected to escape more easily out of low-mass galaxies, and to expel a 
greater fraction of the gas from dwarf galaxies. At high redshifts, most galax- 
ies form in relatively low-mass halos, and the high halo merger rate leads to 
vigorous star formation. Thus, outflows may have had a great impact on the 
earliest generations of galaxies, with conscqncnicx^s that may include metal en- 
richment of the IGM and the disruption of dwarf galaxies. In this subsection 
we present a simple model for the propagation of individual supernova shock 
fronts in the IGM. Wc discuss some implications of this model, but wc defer 
to the following subsection the brunt of the discussion of the cosmological 
consequences of outflows. 

For a galaxy forming in a given halo, the supernova rate is related to the 
star formation rate. In particular, for a Scalo (1998) [315] initial stellar mass 
function, if we assume that a supernova is produced by each M > 8Mq star, 
then on average one supernova explodes for every 126 Mq of star formation, 
expelling an ejecta mass of ~ 3 Mq including ~ 1 Mq of heavy elements. We 
assume that the individual supernovae produce expanding hot bubbles which 
merge into a single overall region delineated by an outwardly moving shock 
front. We assume that most of the baryons in the outflow lie in a thin shell, 
while most of the thermal energy is carried by the hot interior. The total 
ejected mass equals a fraction /gas of the total halo gas which is lifted out 
of the halo by the outflow. This gas mass includes a fraction /eject of the 
mass of the supernova ejecta itself (with /eject < 1 since some metals may be 
deposited in the disk and not ejected). Since at high rcdshift most of the halo 
gas is likely to have cooled onto a disk, we assume that the mass carried by the 
outflow remains constant until the shock front reaches the halo virial radius. 
We assume an average supernova energy of 10^"'^i?5i erg, a fraction /wind of 
which remains in the outflow after it escapes from the disk. The outflow must 
overcome the gravitational potential of the halo, which we assume to have 
a Navarro, Frcnk, & White (1997) [266] density profile [NEW; sec equation 
(88)]. Since the entire shell mass must be lifted out of the halo, we include 
the total shell mass as well as the total injected energy at the outset. This 
assumption is consistent with the fact that the burst of star formation in a 
halo is typically short compared to the total time for which the corresponding 
outflow expands. 

The escape of an outflow from an NFW halo depends on the concentration 
parameter cn of the halo. Simulations by Bullock et al. (2000) [72] indicate 
that the concentration parameter decreases with redshift, and their results 
may be extrapolated to our regime of interest (i.e., to smaller halo masses 
and higher redshifts) by assuming that 



Although we calculate below the dynamics of each outflow in detail, it is also 

useful to estimate which halos can generate large-scale outflows by comparing 
the kinetic energy of the outflow to the potential energy needed to completely 




(139) 



First Light 103 



escape (i.e., to infinite distance) from an NFW halo. We thus find that the 
outflow can escape from its originating halo if the circular velocity is below a 
critical value given by 



T/ onn /^5l/wind('7/0-l) , -1 

Vcrit = 200 W — r km s , (140) 

V /gas S'(cn) 

where the efficiency r/ is the fraction of baryons incorporated in stars, and 

= (l + x)ln(l + a.)-a. " ^'^'^ 

Note that the contribution to /gas of the supernova ejecta itself is 0.024r7/oject) 
so the ejecta mass is usually negligible unless /gas < 1%. Equation (140) can 
also be used to yield the maximum gas fraction /gas which can be ejected from 
halos, as a function of their circular velocity. Although this equation is most 
general, if we assume that the parameters /gas and /wind ai"© independent of 
M and z then we can normalize them based on low-redshift observations. If 
we specify cn ~ 10 (with ^(10) = 6.1) at ^ = 0, then setting E^i = 1 and 
T} = 10% yields the required energy efficiency as a function of the ejected halo 
gas fraction: 

Krit 



/wind — l*^/g 



100 km s"^ 



(142) 



A value of V|;rit ^ 100 km s""'^ is suggested by several theoretical and ob- 
servational arguments which are discussed in the next subsection. However, 
these arguments are not conclusive, and V^rit may differ from this value by 
a large factor, especially at high redshift (where outflows are observationally 
unconstrained at present). Note the degeneracy between /gas and /wind which 
remains even if V^rit is specified. Thus, if 14rit ~ 100 km s~^ then a high 
efficiency /wind ^ 1 is required to eject most of the gas from all halos with 
Vc < Vcvit, but only /wind ~ 10% is required to eject 5-10% of the gas. The 
evolution of the outflow does depend on the value of /wind and not just the 
ratio /wind//gas, sincc the shell accumulates material from the IGM which 
eventually dominates over the initial mass carried by the outflow. 

We solve numerically for the spherical expansion of a galactic outflow, 
elaborating on the basic approach of Tegmark, Silk, & Evrard (1993) [358]. 
We assume that most of the mass m carried along by the outflow lies in a 
thin, dense, relatively cool shell of proper radius R. The interior volume, while 
containing only a fraction /jnt <C 1 of the mass m, carries most of the thermal 
energy in a hot, isothermal plasma of pressure pint and temperature T. We 
assume a uniform exterior gas, at the mean density of the Universe (at each 
redshift), which may be neutral or ionized, and may exert a prcssiire poxt as 
indicated below. We also assume that the dark matter distribution follows the 
NFW profile out to the virial radius, and is at the mean density of the Universe 
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outside the halo virial radius. Note that in reality an overdense distribution 

of gas as well as dark matter may surround each halo due to secondary infall. 
The shell radius R in general evolves as follows: 

(143) 

where the right-hand-side includes forces due to pressure, sweeping up of 
additional mass, gravity, and a cosmological constant, respectively. The shell 
is accelerated by internal pressure and decelerated by external pressure, i.e., 
Sp = Pint — Pext- In the gravitational force, M{R) is the total enclosed mass, 
not including matter in the shell, and is the effective contribution of 
the shell mass in the thin-shell approximation [279]. The interior pressure is 
determined by energy conservation, and evolves according to [358]: 

dpint ^ L _r^PintdR 

dt 2ttR^ R dt ' ^ ' 

where the luminosity L incorporates heating and cooling terms. We include 
in L the supernova luminosity ign (during a brief initial period of energy 
injection), cooling terms Lcoob ionization Lion- f^d dissipation -Ldiss- For sim- 
plicity, we assume ionization equilibrium for the interior plasma, and a pri- 
mordial abundance of hydrogen and helium. We include in icooi all relevant 
atomic cooling processes in hydrogen and helium, i.e., coUisional processes, 
Bremsstrahlung emission, and Compton cooling off the CMB. Compton scat- 
tering is the dominant cooling process for high-redshift outflows. We include 
in Lion only the power required to ionize the incoming hydrogen upstream, at 
the energy cost of 13.6 eV per hydrogen atom. The interaction between the 
expanding shell and the swept-up mass dissipates kinetic energy. The fraction 
fd of this energy which is re- injected into the interior depends on complex pro- 
cesses occurring near the shock front, including turbulence, non-equilibrium 
ionization and cooling, and so (following Tegmark et al. 1993 [358]) we let 

where we set /d = 1 and compare below to the other extreme of fd = 0. 

In an expanding Universe, it is preferable to describe the propagation of 
outflows in terms of comoving coordinates since, e.g., the critical result is the 
maximum comoving size of each outflow, since this size yields directly the 
total IGM mass which is displaced by the outflow and injected with metals. 
Specifically, we apply the following transformation [328, 374]: 

dt = ar'^dt, R = a~^R, p = aPp, p = a?p . (146) 

For flA = 0, Voit (1996) [374] obtained (with the time origin £ = at redshift 
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t = 



QrnHl 



a/1 + fimZl - ^1 + fir, 



(147) 



while for Qm + f^A = 1 there is no simple analytic expression. We set (3 = 
R/fvii, in terms of the virial radius fvir [equation (84)] of the source halo. We 
define ag as the ratio of the shell mass m to ^npb f^jj., where pb = pi,{z = 0) is 
the mean baryon density of the Universe at 2 = 0. More generally, we define 



gTrpf, p-^ 



if /?< 1 
l) I (3^ otherwise. 



(148) 



Here we assumed, as noted above, that the shell mass is constant until the 
halo virial radius is reached, at which point the outflow begins to sweep up 
material from the IGM. We thus derive the following equations: 



(PR 
dp 



as(/3) Pb R 



as(l3)R 



\dt) 



m 



if/3< 1 



^RH^OmS{l3) + lRH^f2bas{l3) otherwise. 



along with 



(149) 



(150) 



In the evolution equation for _R, for /? < 1 we assume for simplicity that the 
baryons are distributed in the same way as the dark matter, since in any case 
the dark matter halo dominates the gravitational force. For /3 > 1, however, 
wc correct (via the last term on the right-hand side) for the presence of mass in 
the shell, since at /? 1 this term may become important. The /? > 1 equation 
also includes the braking force due to the swept-up IGM mass. The enclosed 
mean overdensity for the NFW profile [Eq. (88)] surrounded by matter at the 
mean density is 



m = { 



ln(l+c)-c/(l+c) 



i( 



ln(l+CN/3)-CN/3/(l+CN/3) if ^ <; 1 

otherwise. 



1 



(151) 
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The physics of supernova shells is discussed in Ostriker & McKee (1988) 
[279] along with a number of analytical solutions. The propagation of cos- 
mological blast waves has also been computed by Ostriker & Cowie (1981) 
[278], Bertschinger (1985) [40] and Carr & Ikeuchi (1985) [74]. Voit (1996) 
[374] derived an exact analytic solution to the fluid equations which, although 
of limited validity, is nonetheless useful for understanding roughly how the 
outflow size depends on several of the parameters. The solution requires an 
idealized case of an outflow which at all times expands into a homogeneous 
IGM. Peculiar gravitational forces, and the energy lost in escaping from the 
host halo, are neglected, cooling and ionization losses are also assumed to be 
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negligible, and the external pressure is not included. The dissipated energy 

is assumed to be retained, i.e., is set equal to unity. Under these condi- 
tions, the standard Sedov self-similar solution [324, 325] generalizes to the 
cosmological case as follows [374] : 



where ^ = 2.026 and Eq = Eo/{l + zi)'^ in terms of the initial (i.e., a.tt = t = 
and z = zi) energy Eq. Numerically, the comoving radius is 



In solving the equations described above, we assume that the shock front 
expands into a pre-ionized region which then recombines after a time de- 
termined by the recombination rate. Thus, the external pressure is included 
initially, it is turned off after the pre-ionized region recombines, and it is then 
switched back on at a lower redshift when the Universe is reionized. When 
the ambient IGM is neutral and the pressure is off, the shock loses energy 
to ionization. In practice we find that the external pressure is unimportant 
during the initial expansion, although it is generally important after reioniza- 
tion. Also, at high redshift ionization losses are much smaller than losses due 
to Compton cooling. In the results shown below, we assume an instantaneous 
reionization at 2; = 9. 

Figure 51 shows the results for a starting redshift z = 15, for a halo of 
mass 5.4 x IO^Mq, stellar mass 8.0 x IO^Mq, comoving fyir = 12 kpc, and 
circular velocity = 20 km/s. We show the shell comoving radius in units 
of the virial radius of the source halo (top panel), and the physical peculiar 
velocity of the shock front (bottom panel). Results are shown (solid curve) 
for the standard set of parameters /int = 0.1, fd = 1, /wind = 75%, and 
/gas = 50%. For comparison, we show several cases which adopt the standard 
parameters except for no cooling (dotted curve), no reionization (short-dashed 
curve), /d = (long-dashed curve), or /wind = 15% and /gas = 10% (dot-short 
dashed curve). When reionization is included, the external pressure halts the 
expanding bubble. We freeze the radius at the point of maximum expansion 
(where dR/di = 0), since in reality the shell will at that point begin to spread 
and fill out the interior volume due to small-scale velocities in the IGM. For the 
chosen parameters, the bubble easily escapes from the halo, but when /wind 
and /gas are decreased the accumulated IGM mass slows down the outflow 
more effectively. In all cases the outflow reaches a size of 10 20 times f^ir, i.e., 
100-200 comoving kpc. If all the metals are ejected (i.e., /eject = !)> then this 
translates to an average metallicity in the shell of ~ l-5x 10~^ in units of the 
solar metallicity (which is 2% by mass). The asymptotic size of the outflow 
varies roughly as /^(f^j, as predicted by the simple solution in equation (152), 




(152) 




(153) 
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but the asymptotic size is rather insensitive to /gas (at a fixed /wind) since the 
outflow mass becomes dominated by the swept- up IGM mass once R > 4fvir • 
With the standard parameter values (i.e., those corresponding to the solid 
curve), Figure 51 also shows (dot-long dashed curve) the Voit (1996) [374] 
solution of equation (152). The Voit solution behaves similarly to the no- 
reionization curve at low redshift, although it overestimates the shock radius 
by ~ 30%, and the overestimate is greater compared to the more realistic case 
which does include reionization. 
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Fig. 51. Evolution of a supernova outflow from a, z = 15 halo of circular velocity 
Vc = 20 km/s. Plotted are the shell comoving radius in units of the virial radius 
of the source halo (top panel), and the physical peculiar velocity of the shock front 
(bottom panel). Results are shown for the standard parameters /int ~ 0.1, fd ~ 1, 
/wind ~ 75%, and /gas = 50% (solid curve). Also shown for comparison are the cases 
of no cooling (dotted curve), no reionization (short-dashed curve), /d = (long- 
dashed curve), or /wind ~ 15% and /gas ~ 10% (dot-short dashed curve), as well as 
the simple Voit (1996) [374] solution of equation (152) for the standard parameter 
set (dot-long dashed curve). In cases where the outflow halts, we freeze the radius 
at the point of maximum expansion. 



Figure 52 shows different curves than Figure 51 but on an identical layout. 
A single curve starting at z = 15 (solid curve) is repeated from Figure 51, 
and it is compared here to outflows with the same parameters but starting at 
z = 20 (dotted curve), z = 10 (short-dashed curve), and z = 5 (long-dashed 
curve). A Vc = 20 km/s halo, with a stellar mass equal to 1.5% of the total 
halo mass, is chosen at the three higher redshifts, but at z = 5 a Vc = 42 km/s 
halo is assumed. Because of the suppression of gas infall after reionization, we 
assume that the z = 5 outflow is produced by supernovae from a stellar mass 
equal to only 0.3% of the total halo mass (with a similarly reduced initial shell 
mass), thus leading to a relatively small final shell radius. The main conclusion 



108 Abraham Loeb 



from both figures is the following: In all cases, the outflow undergoes a rapid 
initial expansion over a fractional redshift interval bzj z ~ 0.2, at which point 
the shell has slowed down to ^ 10 km/s from an initial 300 km/s. The rapid 
deceleration is due to the accumulating IGM mass. External pressure from the 
reionized IGM completely halts all high-redshift outflows, and even without 
this effect most outflows would only move at ~ 10 km/s after the brief initial 
expansion. Thus, it may be possible for high-redshift outflows to pollute the 
Lyman alpha forest with metals without affecting the forest hydrodynamically 
at z < 4. While the bulk velocities of these outflows may dissipate quickly, 
the outflows do sweep away the IGM and create empty bubbles. The resulting 
effects on observations of the Lyman alpha forest should be studied in detail 
(some observational signatures of feedback have been suggested recently by 
Theuns, Mo, & Schaye 2000 [362]). 




Fig. 52. Evolution of supernova outflows at different redshifts. The top and bot- 
tom panels are arranged similarly to Figure 51. The z = Vh outflow (solid curve) is 
repeated from Figure 51, and it is compared here to outflows with the same param- 
eters but starting at 2 = 20 (dotted curve), z = 10 (short-dashed curve), and 2 = 5 
(long-dashed curve). A 14 = 20 km/s halo is assumed except for 2 = 5, in which 
case a,Vc~42 km/s halo is assumed to produce the outflow (see text). 



Furlanetto & Loeb (2003) [141] derived the evolution of the characteristic 
scale and filling fraction of supernova-driven bubbles based on a refinement 
of this formalism (see also their 2001 paper for quasar-driven outflows). The 
role of metal-rich outflows in smearing the transition epoch between Pop-Ill 
(metal-free) and Pop II (metal-enriched) stars, was also analysed by Furlan- 
etto & Loeb (2005) [144], who concluded that a double-reionization history in 
which the ionization fraction goes through two (or more) peaks is unlikely. 
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8.2 Effect of Outflows on Dwarf Galaxies and on the IGM 

Galactic outflows represent a complex feedback process which affects the evo- 
lution of cosmic gas through a variety of phenomena. Outflows inject hydro- 
dynamic energy into the interstellar medium of their host galaxy. As shown in 
the previous subsection, even a small fraction of this energy suffices to eject 
most of the gas from a dwarf galaxy, perhaps quenching further star forma- 
tion after the initial burst. At the same time, the enriched gas in outflows 
can mix with the interstellar medium and with the surrounding IGM, allow- 
ing later generations of stars to form more easily because of metal-enhanced 
cooling. On the other hand, the expanding shock waves may also strip gas in 
surrounding galaxies and suppress star formation. 

Dekel & Silk (1986) [104] attempted to explain the different properties of 
diffuse dwarf galaxies in terms of the effect of galactic outflows. They noted 
the observed trends whereby lower-mass dwarf galaxies have a lower surface 
brightness and metallicity, but a higher mass-to- light ratio, than higher mass 
galaxies. They argued that these trends are most naturally explained by sub- 
stantial gas removal from an underlying dark matter potential. Galaxies lying 
in small halos can eject their remaining gas after only a tiny fraction of the 
gas has turned into stars, while larger galaxies require more substantial star 
formation before the resulting outflows can expel the rest of the gas. Assuming 
a wind eflaciency /wind ~ 100%, Dekel & Silk showed that outflows in halos 
below a circular velocity threshold of Vcrit ^ 100 km/s have sufficient energy 
to expel most of the halo gas. Furthermore, cooling is very efficient for the 
characteristic gas temperatures associated with V^rit < 100 km/s halos, but 
it becomes less efficient in more massive halos. As a result, this critical veloc- 
ity is expected to signify a dividing line between bright galaxies and diffuse 
dwarf galaxies. Although these simple considerations may explain a number 
of observed trends, many details are still not conclusively determined. For 
instance, even in galaxies with sufficient energy to expel the gas, it is possible 
that this energy gets deposited in only a small fraction of the gas, leaving the 
rest almost unaffected. 

Since supernova explosions in an inhomogcneous interstellar medium lead 
to complicated hydrodynamics, in principle the best way to determine the 
basic parameters discussed in the previous subsection (/wind. /gas. and /eject) 
is through detailed numerical simulations of individual galaxies. Mac Low & 
Ferrara (1999) [235] simulated a gas disk within a z dark matter halo. 
The disk was assumed to be azimuthally symmetric and initially smooth. They 
represented supernovae by a central source of energy and mass, assuming a 
constant luminosity which is maintained for 50 million years. They found that 
the hot, metal-enriched ejecta can in general escape from the halo much more 
easily than the colder gas within the disk, since the hot gas is ejected in a tube 
perpendicular to the disk without displacing most of the gas in the disk. In 
particular, most of the metals were expelled except for the case with the most 
massive halo considered (with lO^M© in gas) and the lowest luminosity (10^'' 
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erg/s, or a total injection of 2 x 10^^ erg). On the other hand, only a small 
fraction of the total gas mass was ejected except for the least massive halo 
(with IO^Mq in gas), where a luminosity of 10'^^ erg/s or more expelled most 
of the gas. We note that beyond the standard issues of numerical resolution 
and convergence, there are several difficulties in applying these results to 
high-redshift dwarf galaxies. Clumping within the expanding shells or the 
ambient interstellar medium may strongly affect both the cooling and the 
hydrodynamics. Also, the effect of distributing the star formation throughout 
the disk is unclear since in that case several characteristics of the problem 
will change; many small explosions will distribute the same energy over a 
larger gas volume than a single large explosion [as in the Sedov (1959) [324] 
solution; see, e.g., equation (152)], and the geometry will be different as each 
bubble tries to dig its own escape route through the disk. Also, high-redshift 
disks should be denser by orders of magnitude than z = disks, due to the 
higher mean density of the Universe at early times. Thus, further numerical 
simulations of this process are required in order to assess its significance during 
the reionization epoch. 

Some input on these issues also comes from observations. Martin (1999) 
[247] showed that the hottest extended X-ray emission in galaxies is charac- 
terized by a temperature of ^ 10^'^ K. This hot gas, which is lifted out of the 
disk at a rate comparable to the rate at which gas goes into new stars, could 
escape from galaxies with rotation speeds of < 130 km/s. However, these 
results are based on a small sample which includes only the most vigorous 
star-forming local galaxies, and the mass-loss rate depends on assumptions 
about the poorly understood transfer of mass and energy among the various 
phases of the interstellar medium. 

Many authors have attempted to estimate the overall cosmological effects 
of outflows by combining simple models of individual outflows with the forma- 
tion rate of galaxies, obtained via semi-analytic methods [98, 358, 374, 265, 
129, 317] or numerical simulations [148, 149, 80, 9]. The main goal of these cal- 
culations is to explain the characteristic metallicities of different environments 
as a function of redshift. For example, the IGM is observed to be enriched with 
metals at redshifts z < 5. Identification of C IV, Si IV an O VI absorption 
lines which correspond to Lya absorption lines in the spectra of high-redshift 
quasars has revealed that the low-density IGM has been enriched to a metal 
abundance (by mass) of Ziqm ~ lO"^-^^^"-^)^©, where Z© = 0.019 is the 
solar metaUicity [252, 372, 347, 229, 99, 346, 121]. The metal enrichment has 
been clearly identified down to H I column densities of ^ lO"'^^'^ cm~^. The 
detailed comparison of cosmological hydrodynamic simulations with quasar 
absorption spectra has established that the forest of Lya absorption lines is 
caused by the smoothly-fluctuating density of the neutral component of the 
IGM [84, 405, 180]. The simulations show a strong correlation between the 
H I column density and the gas overdensity 6gas [102], implying that metals 
were dispersed into regions with an overdensity as low as ^gas ^ 3 or possibly 
even lower. 
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In general, dwarf galaxies are expected to dominate metal enrichment at 
high-rcdshift for several reasons. As noted above and in the previous sub- 
section, outflows can escape more easily out of the potential wells of dwarfs. 
Also, at high redshift, massive halos are rare and dwarf halos are much more 
common. Finally, as already noted, the Sedov (1959) [324] solution [or equa- 
tion (152)] implies that for a given total energy and expansion time, multiple 
small outflows fill large volumes more effectively than would a smaller number 
of large outflows. Note, however, that the strong cflect of feedback in dwarf 
galaxies may also quench star formation rapidly and reduce the efficiency of 
star formation in dwarfs below that found in more massive galaxies. 

Ccn & Ostrikcr (1999) [80] showed via numerical simulation that metals 
produced by supernovae do not mix uniformly over cosmological volumes. In- 
stead, at each epoch the highest density regions have much higher metallicity 
than the low-density IGM. They noted that early star formation occurs in the 
most overdense regions, which therefore reach a high metallicity (of order a 
tenth of the solar value) by 2: ~ 3, when the IGM metallicity is lower by 1-2 or- 
ders of magnitude. At later times, the formation of high-temperature clusters 
in the highest-density regions suppresses star formation there, while lower- 
density regions continue to increase their metallicity. Note, however, that the 
spatial resolution of the hydrodynamic code of Cen & Ostriker is a few hun- 
dred kpc, and anything occurring on smaller scales is inserted directly via 
simple parametrized models. Scannapieco & Broadhurst (2000) [317] imple- 
mented expanding outflows within a numerical scheme which, while not a full 
gravitational simulation, did include spatial correlations among halos. They 
showed that winds from low-mass galaxies may also strip gas from nearby 
galaxies (sec also Scannapieco, Ferrara, & Broadhurst 2000 [318]), thus sup- 
pressing star formation in a local neighborhood and substantially reducing the 
overall abundance of galaxies in halos below a mass of IO^^Mq. Although 
quasars do not produce metals, they may also affect galaxy formation in their 
vicinity via energetic outflows [116, 15, 339, 263]. 

Gnedin & Ostriker (1997) [148] and Gnedin (1998) [149] identified another 
mixing mechanism which, they argued, may be dominant at high redshift 
(2; > 4). In a collision between two protogalaxies, the gas components collide 
in a shock and the resulting pressure force can eject a few percent of the gas 
out of the merger remnant. This is the merger mechanism, which is based 
on gravity and hydrodynamics rather than direct stellar feedback. Even if 
supernovae inject most of their metals in a local region, larger-scale mixing 
can occur via mergers. Note, however, that Gnedin's (1998) [149] simulation 
assumed a comoving star formation rate at 2; > 5 of ~ ^Mq per year per 
comoving Mpc^, which is 5-10 times larger than the observed rate at redshift 
3 4. Aguirre et al. [9] used outflows implemented in simulations to conclude 
that winds of ~ 300 km/s at 2 < 6 can produce the mean metallicity observed 
at 2 ~ 3 in the Lya forest. In a separate paper Aguirre et al. [10] explored 
another process, where metals in the form of dust grains are driven to large 
distances by radiation pressure, thus producing large-scale mixing without 
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displacing or heating large volumes of IGM gas. The success of this mechanism 
depends on detailed microphysics such as dust grain destruction and the effect 
of magnetic fields. The scenario, though, may be directly testable because it 
leads to significant ejection only of elements which solidify as grains. 

Feedback from galactic outflows encompasses a large variety of processes 
and influences. The large range of scales involved, from stars or quasars em- 
bedded in the interstellar medium up to the enriched IGM on cosmological 
scales, make possible a multitude of different, complementary approaches, 
promising to keep galactic feedback an active field of research. 

9 The Frontier of 21cm Cosmology 
9.1 Mapping Hydrogen Before Reionization 

The small residual fraction of free electrons after cosmological recombination 

coupled the temperature of the cosmic gas to that of the cosmic microwave 
background (CMB) down to a redshift, z ~ 200 [284]. Subsequently, the gas 
temperature dropped adiabatically as Tgas oc (1 + zf' below the CMB tem- 
perature cx (1 + z). The gas heated up again after being exposed to the 
photo-ionizing ultraviolet light emitted by the first stars during the reioniza- 
tion epoch at ^ < 20. Prior to the formation of the first stars, the cosmic 
neutral hydrogen must have resonantly absorbed the CMB flux through its 
spin-flip 21cm transition [131, 323, 367, 404]. The linear density fluctuations 
at that time should have imprinted anisotropies on the CMB sky at an ob- 
served wavelength of A = 21.12[(1 + z)/100] meters. We discuss these early 
21cm fluctuations mainly for pedagogical purposes. Detection of the earliest 
21cm signal will be particularly challenging because the foreground sky bright- 
ness rises as A^ ''' at long wavelengths in addition to the standard \[\ scaling 
of the detector noise temperature for a given integration time and fractional 
bandwidth. The discussion in this section follows Loeb & Zaldarriaga (2004) 



We start by calculating the history of the spin temperature, T,, defined 
through the ratio between the number densities of hydrogen atoms in the 
excited and ground state levels, ni/no = (ffi/ffo) exp {— T^/Tg} , 



where subscripts 1 and correspond to the excited and ground state levels of 
the 21cm transition, [gi/go) = 3 is the ratio of the spin degeneracy factors of 
the levels, tih = {no + n\) oc (1 -|- z)'^ is the total hydrogen density, and = 
0.068K is the temperature corresponding to the energy difference between the 
levels. The time evolution of the density of atoms in the ground state is given 
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Fig. 53. The 21cm transition of hydrogen. The higher energy level the spin of 
the electron (e-) is aligned with that of the proton (p+). A spin flip results in the 
emission of a photon with a wavelength of 21cm (or a frequency of 1420MHz). 



where a{t) = (1 + z)^^ is the cosmic scale factor, yl's and B's are the Einstein 
rate coefficients, C's are the collisional rate coefficients, and is the black- 
body intensity in the Rayleigh- Jeans tail of the CMB, namely = 2kT^/}? 
with A = 21 cm [306]. Here a dot denotes a time-derivative. The 0^1 transi- 
tion rates can be related to the 1 — > transition rates by the requirement that 
in thermal equilibrium with Tg = = Tgas, the right-hand-side of Eq. (155) 
should vanish with the collisional terms balancing each other separately from 
the radiative terms. The Einstein coefficients are Aiq — 2.85 x 10^^^ s~^, 
Bw = {X^/2hc)Aio and B^i = (51/30)^10 [131, 306]. The collisional de- 
excitation rates can be written as Cio — |'«(1 — 0)nH, where k(1 — 0) is 
tabulated as a function of Tgas [11, 406]. 

Equation (155) can be simplified to the form. 




(155) 



dz 



+ (l-r)(Cio + Aio + Sio/.)], 



(156) 
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where T = no/nn, H w Hq^/IT^{1 + z)'^/^ is the Hubble parameter at high 
redshifts (with a present-day value of ifo), and flm is the density parameter of 
matter. The upper panel of Fig. 54 shows the results of integrating Eq. (156). 
Both the spin temperature and the kinetic temperature of the gas track the 
CMB temperature down to z ~ 200. Collisions are efficient at coupling T, 
and Tgas down to z ~ 70 and so the spin temperature follows the kinetic tem- 
perature around that redshift. At much lower redshifts, the Hubble expansion 
makes the collision rate subdominant relative the radiative coupling rate to 
the CMB, and so Tg tracks again. Consequently, there is a redshift window 
between 30 < z < 200, during which the cosmic hydrogen absorbs the CMB 
flux at its resonant 21cm transition. Coincidentally, this redshift interval pre- 
cedes the appearance of collapsed objects [23] and so its signatures are not 
contaminated by nonlinear density structures or by radiative or hydrodynamic 
feedback effects from stars and quasars, as is the case at lower redshifts [404]. 




1 10 100 1000 



(1+z) 

Fig. 54. Upper panel: Evolution of the gas, CMB and spin temperatures with red- 
shift [4]. Lower panel: dTt/dSn as function of redshift. The separate contributions 
from fluctuations in the density and the spin temperature are depicted. We also 
show dTfj/dfeo oc dT^/dSH x 5h, with an arbitrary normalization. 



During the period when the spin temperature is smaller than the CMB 
temperature, neutral hydrogen atoms absorb CMB photons. The resonant 
21cm absorption reduces the brightness temperature of the CMB by, 

n = T{T,-T.,)/{\ + z), (157) 

where the optical depth for resonant 21cm absorption is. 



32nkTsH{z) ' 



(158) 
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Small inhomogeneities in the hydrogen density Sh = (wh — ^h)/^h result 
in fluctuations of the 21cm absorption through two separate effects. An excess 
of neutral hydrogen directly increases the optical depth and also alters the 
evolution of the spin temperature. For now, we ignore the additional effects 
of pecuhar velocities (Bharadwaj & Ah 2004 [41]: Barkana & Locb 2004 [27]) 
as well as fluctuations in the gas kinetic temperature due to the adiabatic 
compression (rarefaction) in overdense (underdense) regions [29] . Under these 
approximations, we can write an equation for the resulting evolution of T 
fluctuations, 

^ = [H{1 + z)]-' {[Cio + Coi + (Boi + B,o)Iu]5r 

az 

+ [Coir-Cw{i-r)]5H], (159) 

leading to spin temperature fluctuations, 

5Ts 1 5T 



T, ln[3r/(l - T)] r(l - Y) ' 



(160) 



The resulting brightness temperature fluctuations can be related to the deriva- 
tive, 

The spin temperature fluctuations STg/Ts are proportional to the density 
fluctuations and so we define, 

dTb _^ , T^fb 5Ts , . 

aOH {rs-Tj)TsdH 

through STb = {dTb/dSH)5H- We ignore fluctuations in dj due to fluctuations 
in Tgas which are very small [11]. Figure 54 shows dTb/ddn as a function of 
redshift, including the two contributions to dTb/dSn, one originating directly 
from density fluctuations and the second from the associated changes in the 
spin temperature [323]. Both contributions have the same sign, because an 
increase in density raises the collision rate and lowers the spin temperature 
and so it allows Tg to better track Tgas- Since 5h grows with time as 6h oc a, 
the signal peaks at 2 ~ 50, a slightly lower redshift than the peak of dTb/dSH- 
Next we calculate the angular power spectrum of the brightness tempera- 
ture on the sky, resulting from density perturbations with a power spectrum 

(Jff (ki)(5ff (k2)) = (27r)35^(ki + k2)P5{h). (163) 

where ^//(k) is the Fourier tansform of the hydrogen density field, k is the 
comoving wavevector, and (• • • ) denotes an ensemble average (following the 
formalism described in [404]). The 21cm brightness temperature observed at 
a frequency v corresponding to a distance r along the line of sight, is given by 
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(164) 



where n denotes the direction of observation, W^{r) is a narrow function of 
r that peaks at the distance corresponding to v. The details of this function 
depend on the characteristics of the experiment. The brightness fluctuations 
in Eq. 164 can be expanded in spherical harmonics with expansion coeffi- 
cients aim{v)- The angular power spectrum of map Ciiv) = {\aim{v)\'^) can 
be expressed in terms of the 3D power spectrum of fluctuations in the density 
ft(fc), 



(2^ 



Ps{k)af{k,iy) 



ai(k, v) 



drWr„ir)-Tr-ir)jiikr). 

dOH 



(165) 



Our calculation ignores inhomogeneities in the hydrogen ionization fraction, 
since they freeze at the earlier recombination epoch (z ~ 10'^) and so their 
amplitude is more than an order of magnitude smaller than 5// at z < 100. 
The gravitational potential perturbations induce a redshift distortion effect 
that is of order ^ (H/ck)'^ smaller than 6h for the high~Z modes of interest 
here. 
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Fig. 55. Angular power spectrum of 21cm anisotropies on the sky at various red- 
shifts. From top to bottom, z = 55, 40, 80, 30, 120, 25, 170. 



Figure 55 shows the angular power spectrum at various redshifts. The sig- 
nal peaks around z ^ 50 but maintains a substantial amplitude over the full 
range of 30 < z < 100. The ability to probe the small scale power of density 
fluctuations is only limited by the Jeans scale, below which the dark matter 
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inhomogeneities are washed out by the finite pressure of the gas. Interestingly, 
the cosmological Jeans mass roaches its minimum value, ^ 3 x 10^ Mq, within 
the redshift interval of interest here which corresponds to modes of angular 
scale ~ arcsecond on the sky. During the epoch of reionization, photoionization 
heating raises the Jeans mass by several orders of magnitude and broadens 
spectral features, thus limiting the ability of other probes of the intergalactic 
medium, such as the Lya forest, from accessing the same very low mass scales. 
The 21cm tomography has the additional advantage of probing the majority 
of the cosmic gas, instead of the trace amount (~ 10~^) of neutral hydrogen 
probed by the Lya forest after reionization. Similarly to the primary CMB 
anisotropics, the 21cm signal is simply shaped by gravity, adiabatic cosmic ex- 
pansion, and well-known atomic physics, and is not contaminated by complex 
astrophysical processes that affect the intergalactic medium at 2; < 30. 

Characterizing the initial fluctuations is one of the primary goals of ob- 
servational cosmology, as it offers a window into the physics of the very early 
Universe, namely the epoch of inflation during which the fluctuations are be- 
lieved to have been produced. In most models of inflation, the evolution of the 
Hubble parameter during inflation leads to departures from a scale-invariant 
spectrum that are of order l/A^efoid with A/efoid ~ 60 being the number of 
enfolds between the time when the scale of our horizon was of order the hori- 
zon during inflation and the end of inflation [218]. Hints that the standard 
ylCDM model may have too much power on galactic scales have inspired sev- 
eral proposals for suppressing the power on small scales. Examples include 
the possibility that the dark matter is warm and it decoupled while being 
relativistic so that its free streaming erased small-scale power [48], or direct 
modifications of inflation that produce a cut-off in the power on small scales 
[192]. An unavoidable collisionless component of the cosmic mass budget be- 
yond CDM, is provided by massive neutrinos (see [198] for a review). Particle 
physics experiments established the mass splittings among different species 
which translate into a lower limit on the fraction of the dark matter accounted 
for by neutrinos of fi, > 0.3%, while current constraints based on galaxies as 
tracers of the small scale power imply < 12% [360]. 

Figure 56 shows the 21cm power spectrum for various models that differ 
in their level of small scale power. It is clear that a precise measurement of 
the 21cm power spectrum will dramatically improve current constraints on 
alternatives to the standard ylCDM spectrum. 

The 21cm signal contains a wealth of information about the initial fluc- 
tuations. A full sky map at a single photon frequency measured up to /max. 
can probe the power spectrum up to /cmax ~ (^max/10'*)Mpc~^. Such a map 
contains l^^y. independent samples. By shifting the photon frequency, one 
may obtain many independent measurements of the power. When measur- 
ing a mode /, which corresponds to a wavenumber k ~ l/r, two maps at 
different photon frequencies will be independent if they are separated in 
radial distance by l/k. Thus, an experiment that covers a spatial range 
Ar can probe a total of kAr ~ lAr/r independent maps. An experiment 
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Fig. 56. Upper panel: Power spectrum of 21cm anisotropies at 2 = 55 for a /ICDM 
scale-invariant power spectrum, a model with n — 0.98, a model with n = 0.98 and 
Or = In P/d In fc^) — —0.07, a model of warm dark matter particles with a mass 
of 1 keV, and a model in which = 10% of the matter density is in three species 
of massive neutrinos with a mass of 0.4 eV each. Lower panel: Ratios between the 
different power spectra and the scale-invariant spectrum. 



that detects the 21cin signal over a range Ai/ centered on a frequency 
is sensitive to Ar/r ^ 0.5(Z\i//i/)(l -I- z)~^^'^, and so it measures a total of 
^"210111 ~ 3 X 10i6(;,„ax/10^)3(Z\i//i/)(z/100)-i/2 independent samples. 

This detection capability cannot be reproduced even remotely by other 
techniques. For example, the primary CMB anisotropies are damped on small 
scales (through the so-called Silk damping), and probe only modes with I < 
3000 (fc < 0.2 Mpc"""^). The total number of modes available in the full sky 
is A'^cmb — "^Imax 2 X 10''(Zmax/3000)^, including both temperature and 
polarization information. 

The sensitivity of an experiment depends strongly on its particular design, 
involving the number and distribution of the antennae for an interferometer. 
Crudely speaking, the uncertainty in the measurement of [1(1 + l)Ci/27r]^/^ is 
dominated by noise, N^,, which is controlled by the sky brightness at the 
observed frequency v [404], 

\5x IQSjy sr-V V 35 / V V /cover/ 

\ to J V50MHzy \ J 

where Zmin is the minimum observable / as determined by the field of view 
of the instruments, ^max is the maximum observable I as determined by the 
maximum separation of the antennae, /cover is the fraction of the array area 
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thats is covered by telescopes, is the observation time and Av is the fre- 
quency range over which the signal can be detected. Note that the assumed 
sky temperature of 0.7 x lO^K at v = 50MHz (corresponding to z ~ 30) is 
more than six orders of magnitude larger than the signal. We have already 
included the fact that several independent maps can be produced by varying 
the observed frequency. The numbers adopted above are appropriate for the 
inner core of the LOFAR array {http://www.lofar.org), planned for initial op- 
eration in 2006. The predicted signal is ~ ImK, and so a year of integration 
or an increase in the covering fraction are required to observe it with LOFAR. 
Other experiments whose goal is to detect 21cm fluctuations from the sub- 
sequent epoch of reionization at z ^ 6 — 12 (when ionized bubbles exist and 
the fluctuations are larger) include the Mileura Wide-Field Array (MWA; 
http://web.haystack.mit.edu/arrays/MWA/), the Primeval Structure Tele- 
scope {PAST; http://arxiv.org/abs/astro-ph/0502029) , and in the more dis- 
tant future the Square Kilometer Array {SKA; http://www.skatelescope.org). 
The main challenge in detecting the predicted signal from higher redshifts 
involves its appearance at low frequencies where the sky noise is high. Pro- 
posed space-based instruments [194] avoid the terrestrial radio noise and the 
increasing atmospheric opacity at < 20 MHz (corresponding to z > 70). 




Fig. 57. Prototype of the tile design for the Mileura Wide-Field Array (MWA) in 
western Australia, aimed at detecting redshifted 21cm from the epoch of reionization. 
Each 4mx4m tile contains 16 dipole antennas operating in the frequency range of 
80-300MHz. Altogether the initial phase of MWA (the so-called "Low-Frequency 
Demostrator" ) will include 500 antenna tiles with a total collecting area of 8000 m^ 
at 150MHz, scattered across a 1.5 km region and providing an angular resolution of 
a few arcminutes. 

The 21cm absorption is replaced by 21cm emission from neutral hydrogen 
as soon as the intergalactic medium is heated above the CMB temperature 
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by X-ray sources during the epoch of reionization [88]. This occurs long be- 
fore rcionization since the required heating requires only a modest amount 
of energy, ~ 10~^ + z)/30], which is three orders of magnitude smaller 

than the amount necessary to ionize the Universe. As demonstrated by Chen 
& Miralda-Escude (2004) [88] , heating due the recoil of atoms as they absorb 
Lya photons [237] is not effective; the Lya color temperature reaches equi- 
librium with the gas kinetic temperature and suppresses subsequent heating 
before the level of heating becomes substantial. Once most of the cosmic hy- 
drogen is reionized at Zreion, the 21cm signal is diminished. The optical depth 
for free-free absorption after reionization, ~ 0.1[(1 -|- 2;reion)/20]^/^, modifies 
only slightly the expected 21cm anisotropics. Gravitational lensing should 
modify the power spectrum [287] at high I, but can be separated as in stan- 
dard CMB studies (see [326] and references therein). The 21cm signal should 
be simpler to clean as it includes the same lensing foreground in independent 
maps obtained at different frequencies. 

The large number of independent modes probed by the 21cm signal would 

— 1/2 

provide a measure of non-Gaussian deviations to a level of ~ -^210111 ' consti- 
tuting a test of the inflationary origin of the primordial inhomogeneities which 
are expected to possess deviations > 10~^ [245]. 

9.2 The Characteristic Observed Size of Ionized Bubbles at the 

End of Reionization 

The first galaxies to appear in the Universe at redshifts z > 20 created ion- 
ized bubbles in the intergalactic medium (IGM) of neutral hydrogen (HI ) 
left over from the Big-Bang. It is thought that the ionized bubbles grew with 
time, surrounded clusters of dwarf galaxies [6 7, 143] and eventually overlapped 
quickly throughout the Universe over a narrow redshift interval near 2; ~ 6. 
This event signaled the end of the reionization epoch when the Universe was 
a billion years old. Measuring the unknown size distribution of the bubbles 
at their final overlap phase is a focus of forthcoming observational programs 
aimed at highly redshifted 21cm emission from atomic hydrogen. In this sub- 
section we follow Wyithe & Loeb (2004) [399] and show that the combined 
constraints of cosmic variance and causality imply an observed bubble size at 
the end of the overlap epoch of ~ 10 physical Mpc, and a scatter in the ob- 
served redshift of overlap along difi'erent lines-of-sight of ~ 0.15. This scatter 
is consistent with observational constraints from recent spectroscopic data on 
the farthest known quasars. This result implies that future radio experiments 
should be tuned to a characteristic angular scale of ~ 0.5° and have a mini- 
mum frequency band-width of ~ 8 MHz for an optimal detection of 21cm flux 
fluctuations near the end of reionization. 

During the reionization epoch, the characteristic bubble size (defined here 
as the spherically averaged mean radius of the H II regions that contain most 
of the ionized volume[143]) increased with time as smaller bubbles combined 
until their overlap completed and the diffuse IGM was reionized. However the 
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Fig. 58. Schematic sketch of the evolution of the kinetic temperature (Tk) and 
spin temperature (Ts) of cosmic hydrogen. Following cosmological recombination at 
z ^ 10^, the gas temperature (orange curve) tracks the CMB temperature (blue 
line; Tj cx {1 + z)) down to z ~ 200 and then declines below it (T^ oc {1 + z)^) 
until the first X-ray sources (accreting black holes or exploding supernovae) heat it 
up well above the CMB temperature. The spin temperature of the 21cm transition 
(red curve) interpolates between the gas and CMB temperatures. Initially it tracks 
the gas temperature through coUisional coupling; then it tracks the CMB through 
radiative coupling; and eventually it tracks the gas temperature once again after 
the production of a cosmic background of UV photons between the Lya and the 
Lyman-limit frequencies that redshift or cascade into the Lya resonance (through 
the Wouthuysen-Field effect [Wouthuysen 1952 [388]; Field 1959 [131]]). Parts of the 
curve are exaggerated for pedagogic purposes. The exact shape depends on astro- 
physical details about the first galaxies, such as their production of X-ray binaries, 
supernovae, nuclear accreting black holes, and their generation of relativistic elec- 
trons in collisionless shocks which produce UV and X-ray photons through inverse- 
Compton scattering of CMB photons. 

largest size of isolated bubbles (fully surrounded by H I boundaries) that can 
be observed is finite, because of the combined phenomena of cosmic variance 
and causality. Figure 61 presents a schematic illustration of the geometry. 
There is a surface on the sky corresponding to the time along different lines- 
of-sight when the diffuse (uncoUapsed) IGM was most recently neutral. We 
refer to it as the Surface of Bubble Overlap (SBO). There are two competing 
sources for fluctuations in the SBO, each of which is dependent on the char- 
acteristic size, i?SBO; of the ionized regions just before the final overlap. First, 
the finite speed of light implies that 21cm photons observed from different 
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21cm Tomography of Ionized Bubbles During Reionization is tike 



Slicing Swiss Cheese 




Observed wavelength fPdislance 
21cm X (1 + z) 



Fig. 59. 21cm imaging of ionized bubbles during tlie epoch of reionization is anal- 
ogous to slicing Swiss cheese. The technique of slicing at intervals separated by the 
typical dimension of a bubble is optimal for revealing different pattens in each slice. 

points along the curved boundary of an H II region must have been emitted 
at different times during the history of the Universe. Second, bubbles on a 
comoving scale R achieve reionization over a spread of redshifts due to cosmic 
variance in the initial conditions of the density field smoothed on that scale. 
The characteristic scale of H II bubbles grows with time, leading to a decline 
in the spread of their formation redshifts [67] as the cosmic variance is averaged 
over an increasing spatial volume. However the 21cm light-travel time across 
a bubble rises concurrently. Suppose a signal 21cm photon which encodes the 
presence of neutral gas, is emitted from the far edge of the ionizing bubble. 
If the adjacent region along the line-of-sight has not become ionized by the 
time this photon reaches the near side of the bubble, then the photon will en- 
counter diffuse neutral gas. Other photons emitted at this lower redshift will 
therefore also encode the presence of diffuse neutral gas, implying that the 
first photon was emitted prior to overlap, and not from the SBO. Hence the 
largest observable scale of H II regions when their overlap completes, corre- 
sponds to the first epoch at which the light crossing time becomes larger than 
the spread in formation times of ionized regions. Only then will the signal 
photon leaving the far side of the HII region have the lowest redshift of any 
signal photon along that line-of-sight. 

The observed spectra of some quasars beyond z 6.1 show a Gunn- 
Peterson trough[163, 127] (Fan et al. 2005 [128]), a blank spectral region at 
wavelengths shorter than Lya at the quasar redshift, implying the presence 
of HI in the diffuse IGM. The detection of Gunn-Peterson troughs indicates 
a rapid change[126, 288, 381] in the neutral content of the IGM at z ~ 6, 
and hence a rapid change in the intensity of the background ionizing flux. 
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Fig. 60. Spectra of 19 quasars with redshifts 5.74 < z < 6.42 from the Sloan Digital 
Sky Survey [128]. For some of the highest-redshift quasars, the spectrum shows no 
transmitted flux shortward of the Lya wavelength at the quasar redshift (the so- 
called "Gunn- Peterson trough"), indicating a non-negligible neutral fraction in the 
IGM (see the analysis of Fan et al. [128] for details). 



This rapid change implies that overlap, and hence the reionization epoch, 
concluded near z ~ 6. The most promising observational probe[404, 259] of the 
reionization epoch is redshifted 21cm emission from intergalactic HI . Future 
observations using low frequency radio arrays (e.g. LOFAR, MWA, and PAST) 
will allow a direct determination of the topology and duration of the phase of 
bubble overlap. In this section we determine the expected angular scale and 
redshift width of the 21cm fluctuations at the SBO theoretically, and show 
that this determination is consistent with current observational constraints. 

We start by quantifying the constraints of causality and cosmic variance. 
First suppose we have an H II region with a physical radius R/{1 + (z)). For 
a 21cm photon, the light crossing time of this radius is 



(A 



,2\l/2 _ 



R 



(167) 



where at the high-redshifts of interest (dz/dt) = —{Ho\/ f2,„)(l-|-z)^/^. Here, c 
is the speed of light, Hq is the present-day Hubble constant, Hm is the present 
day matter density parameter, and (z) is the mean redshift of the SBO. Note 
that when discussing this crossing time, we are referring to photons used to 
probe the ionized bubble (e.g. at 21cm), rather than photons involved in the 
dynamics of the bubble evolution. 

Second, overlap would have occurred at different times in different regions 
of the IGM due to the cosmic scatter in the process of structure formation 
within finite spatial volumes [67]. Reionization should be completed within 
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a region of comoving radius R when the fraction of mass incorporated into 
collapsed objects in this region attains a certain critical value, corresponding 
to a threshold number of ionizing photons emitted per baryon. The ionization 
state of a region is governed by the enclosed ionizing luminosity, by its over- 
density, and by dense pockets of neutral gas that arc self shielding to ionizing 
radiation. There is an offset [67] 6z between the redshift when a region of 
mean over-density achieves this critical collapsed fraction, and the redshift 
z when the Universe achieves the same collapsed fraction on average. This 
offset may be computed[67] from the expression for the collapsed fraction[52] 
Fcoi within a region of over-density on a comoving scale R, 



Fco\{M^in) = erfc 



(1 + z) Sc{z) 



1 



(168) 

where 5c{z) oc (1 + 2;) is the collapse threshold for an over-density at a redshift 
z: (7r and <7i^^^.^_^ are the variances in the power-spectrum linearly extrapolated 
to z = on comoving scales corresponding to the region of interest and to 
the minimum galaxy mass Mmin, respectively. The offset in the ionization 
redshift of a region depends on its linear over-density, (5r. As a result, the 
distribution of offsets, and therefore the scatter in the SBO may be obtained 
directly from the power spectrum of primordial inhomogeneities. As can be 
seen from equation (168), larger regions have a smaller scatter due to their 
smaller cosmic variance. 

Note that equation (168) is independent of the critical value of the col- 
lapsed fraction required for reionization. Moreover, our mimerical constraints 
are very weakly dependent on the minimum galaxy mass, which we choose 
to have a virial temperature of lO^K corresponding to the cooling threshold 
of primordial atomic gas. The growth of an H II bubble around a cluster 
of sources requires that the mean-free-path of ionizing photons be of order 
the bubble radius or larger. Since ionizing photons can be absorbed by dense 
pockets of neutral gas inside the H II region, the necessary increase in the 
mean-free-path with time implies that the critical collapsed fraction required 
to ionize a region of size R increases as well. This larger collapsed fraction af- 
fects the redshift at which the region becomes ionized, but not the scatter in 
redshifts from place to place which is the focus of this sub-section. Our results 
are therefore independent of assumptions about unknown quantities such as 
the star formation efficiency and the escape fraction of ionizing photons from 
galaxies, as well as unknown processes of feedback in galaxies and clumping 
of the IGM. 

Figure 62 displays the above two fundamental constraints. The causality 
constraint (Eq. 167) is shown as the blue line, giving a longer crossing time 
for a larger bubble size. This contrasts with the constraint of cosmic variance 
(Eq. 168), indicated by the red line, which shows how the scatter in formation 
times decreases with increasing bubble size. The scatter in the SBO redshift 
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and the corresponding fluctuation scale of the SBO are given by the intersec- 
tion of these curves. We find that the thickness of the SBO is (Z\z^)^/^ ^ 0.13, 
and that the bubbles which form the SBO have a characteristic comoving size 
of ~ 60Mpc (equivalent to 8.6 physical Mpc). At 2; ~ 6 this size corresponds 
to angular scales of 0sbo ~ 0.4 degrees on the sky. 

A scatter of ~ 0.15 in the SBO is somewhat larger than the value ex- 
tracted from existing numerical simulations[152, 402]. The difference is most 
likely due to the limited size of the simulated volumes; while the simulations 
appropriately describe the reionization process within limited regions of the 
Universe, they are not sufSciently large to describe the global properties of the 
overlap phase [67]. The scales over which cosmological radiative transfer has 
been simulated are smaller than the characteristic extent of the SBO, which 
we find to be i?sBO ~ 70 comoving Mpc. 

We can constrain the scatter in the SBO rcdshift obscrvationally using 
the spectra of the highest redshift quasars. Since only a trace amount of neu- 
tral hydrogen is needed to absorb Lya photons, the time where the IGM 
becomes Lya transparent need not coincide with bubble overlap. Following 
overlap the IGM was exposed to ionizing sources in all directions and the 
ionizing intensity rose rapidly. After some time the ionizing background flux 
was sufficiently high that the HI fraction fell to a level at which the IGM 
allowed transmission of resonant Lya photons. This is shown schematically 
in Figure 61. The lower wavelength limit of the Gunn-Peterson trough corre- 
sponds to the Lya wavelength at the redshift when the IGM started to allow 
transmission of Lya photons along that particular line-of-sight. In addition 
to the SBO we therefore also define the Surface of Lya Transmission (here- 
after SLT) as the redshift along different lines-of-sight when the diffuse IGM 
became transparent to Lya photons. 

The scatter in the SLT redshift is an observable which we would like to 
compare with the scatter in the SBO redshift. The variance of the density 
field on large scales results in the biased clustering of sources [6 7]. H II regions 
grow in size around these clusters of sources. In order for the ionizing photons 
produced by a cluster to advance the walls of the ionized bubble around 
it, the mean-free-path of these photons must be of order the bubble size or 
larger. After bubble overlap, the ionizing intensity at any point grows until 
the ionizing photons have time to travel across the scale of the new mean-free- 
path, which represents the horizon out to which ionizing sources are visible. 
Since the mean- free-path is larger than Rsbo, the ionizing intensity at the 
SLT averages the cosmic scatter over a larger volume than at the SBO. This 
constraint implies that the cosmic variance in the SLT redshift must be smaller 
than the scatter in the SBO redshift. However, it is possible that opacity from 
small-scale structure contributes additional scatter to the SLT redshift. 

If cosmic variance dominates the observed scatter in the SLT redshift, 
then based on the spectra of the three 2; > 6.1 quasars[127, 381] we would 
expect the scatter in the SBO redshift to satisfy (Z\z^)^(^g > 0.05. In addition, 
analysis of the proximity effect for the size of the H II regions around the two 
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highest redshift quasars [396, 251] implies a neutral fraction that is of order 
unity (i.e. prc-ovcrlap) at z ^ 6.2 — 6.3, while the transmission of Lya photons 
at z < 6 implies that overlap must have completed by that time. This restricts 
the scatter in the SBO to be {Az^)^Q. < 0.25. The constraints on values for 
the scatter in the SBO redshift are shaded gray in Figure 62. It is reassuring 
tliat the theoretical prediction for the SBO scatter of {Az'^)^^Q, ^ 0.15, with a 
characteristic scale of ^ 70 comoving Mpc, is bounded by these constraints. 

The possible presence of a significantly neutral IGM just beyond the red- 
shift of overlap[396, 251] is encouraging for upcoming 21cm studies of the 
reionization epoch as it results in emission near an observed frequency of 200 
MHz where tlie signal is most readily detectable. Future observations of red- 
shifted 21cm line emission at 6 < z < 6.5 with instruments such as LOFAR, 
MWA, and PAST, will be able to map the three-dimensional distribution of HI 
at the end of reionization. The intergalactic H II regions will imprint a 'knee' 
in the power-spectrum of the 21cm anisotropics on a characteristic angular 
scale corresponding to a typical isolated H II region[404]. Our results sug- 
gest that this characteristic angular scale is large at the end of reionization, 
^SBO ~ 0.5 degrees, motivating the construction of compact low frequency 
arrays. An SBO thickness of {Az'^^^ 0.15 suggests a minimum frequency 
band-width of ~ 8 MHz for experiments aiming to detect anisotropics in 21cm 
emission just prior to overlap. These results will help guide the design of the 
next generation of low-frequency radio observatories in the search for 21cm 
emission at the end of the reionization epoch. 

The full size distribution of ionized bubbles has to be calculated from a 
numerical cosmological simulation that includes gas dynamics and radiative 
transfer. The simulation box needs to be sufficiently large for it to sample an 
unbiased volume of the Universe with little cosmic variance, but at the same 
time one must resolve the scale of individual dwarf galaxies which provide (as 
well as consume) ionizing photons (see discussion at the last section of this 
review). Until a reliable simulation of this magnitude exists, one must adopt 
an approximate analytic approach to estimate the bubble size distribution. 
Below we describe an example for such a method, developed by Furlanetto, 
Zaldarriaga, & Hernquist (2004) [143]. 

The criterion for a region to be ionized is that galaxies inside of it produce 
a sufficient number of ionizing photons per baryon. This condition can be 
translated to the requirement that the collapsed fraction of mass in halos above 
some threshold mass Mj^in will exceed some threshold, namely Fco\ > C~^- 
The minimum halo mass most likely corresponds to a virial temperature of 
lO'^K relating to the threshold for atomic cooling (assuming that molecular 
hydrogen cooling is suppressed by the UV background in the Lyman- Werner 
band). We would like to find the largest region around every point that satisfies 
the above condition on the collapse fraction and then calculate the abundance 
of ionized regions of this size. Different regions have different values of Fcoi 
because their mean density is different. In the extended Press-Schechter model 
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Fig. 61. The distances to the observed Surface of Bubble Overlap (SBO) and Surface 
of Lya Transmission (SLT) fluctuate on the sky. The SBO corresponds to the first 
region of diffuse neutral IGM observed along a random line-of-sight. It fluctuates 
across a shell with a minimum width dictated by the condition that the light crossing 
time across the characteristic radius -Rsbo of ionized bubbles equals the cosmic 
scatter in their formation times. Thus, causality and cosmic variance determine 
the characteristic scale of bubbles at the completion of bubble overlap. After some 
time delay the IGM becomes transparent to Lya photons, resulting in a second 
surface, the SLT. The upper panel illustrates how the lines-of-sight towards two 
quasars (Ql in red and Q2 in blue) intersect the SLT with a redshift difference 
5z. The resulting variation in the observed spectrum of the two quasars is shown 
in the lower panel. Observationally, the ensemble of redshifts down to which the 
Gunn-Peterson troughs are seen in the spectra of « > 6.1 quasars is drawn from the 
probability distribution dP/dzsur for the redshift at which the IGM started to allow 
Lya transmission along random lines-of-sight. The observed values of zslt show a 
small scatter[127] in the SLT redshift around an average value of (zslt) ~ 5.95. 
Some regions of the IGM may have also become transparent to Lya photons prior to 
overlap, resulting in windows of transmission inside the Gunn-Peterson trough (one 
such region may have been seen[381] in SDSS J1148-I-5251). In the existing examples, 
the portions of the Universe probed by the lower end of the Gunn-Peterson trough 
are located several hundred comoving Mpc away from the background quasar, and 
are therefore not correlated with the quasar host galaxy. The distribution dP/dzsLT 
is also independent of the redshift distribution of the quasars. Moreover, lines-of- 
sight to these quasars are not causally connected at z ~ 6 and may be considered 
independent. 



128 



Abraham Loeb 




o 
m 
to 



0.1 



1 



<Az^> 



,2^1/2 



Fig. 62. Constraints on the scatter in the SBO redshift and the characteristic size of 
isolated bubbles at the final overlap stage, iisBO (see Fig. 1). The characteristic size 
of H 11 regions grows with time. The SBO is observed for the bubble scale at which 
the light crossing time (blue line) first becomes smaller than the cosmic scatter in 
bubble formation times (red line). At z 6, the implied scale i?sBO ~ 60 comoving 
Mpc (or 8.6 physical Mpc), corresponds to a characteristic angular radius of 
^SLT ~ 0.4 degrees on the sky. After bubble overlap, the ionizing intensity grows to a 
level at which the IGM becomes transparent to Lya photons. The collapsed fraction 
required for Lya transmission within a region of a certain size will be larger than 
required for its ionization. However, the scatter in equation (168) is not sensitive to 
the collapsed fraction, and so may be used for both the SBO and SLT. The scatter 
in the SLT is smaller than the cosmic scatter in the structure formation time on 
the scale of the mean-free-path for ionizing photons. This mean-free-path must be 
longer than TisBO ~ 60Mpc, an inference which is supported by analysis of the Lya 
forest at z ~ 4 where the mean-free-path is estimated [257] to be ~ 120 comoving 
Mpc at the Lyman limit (and longer at higher frequencies). If it is dominated by 
cosmic variance, then the scatter in the SLT redshift provides a lower limit to the 
SBO scatter. The three known quasars at z > 6.1 have Lya transmission redshifts 
of[381, 127] zsLT = 5.9, 5.95 and 5.98, implying that the scatter in the SBO must 
be > 0.05 (this scatter may become better known from follow-up spectroscopy of 
Gamma Ray Burst afterglows at z > 6 that might be discovered by the SWIFT 
satellite [26, 61]). The observed scatter in the SLT redshift is somewhat smaller 
than the predicted SBO scatter, confirming the expectation that cosmic variance is 
smaller at the SLT. The scatter in the SBO redshift must also be < 0.25 because the 
lines-of-sight to the two highest redshift quasars have a redshift of Lya transparency 
at z ~ 6, but a neutral fraction that is known from the proximity effect[396] to be 
substantial at z > 6.2 — 6.3. The excluded regions of scatter for the SBO are shown 
in gray. 
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(Bond et al. 1991 [52]; Lacey & Cole 1993 [212]), the collapse fraction in a 
region of mean overdensity 6m is 

i^coi = erfc ( , ~ I . (169) 

where (t'^{M, z) is the variance of density fluctuations on mass scale M, 0-^;^^ = 

a'^{M„^in, z), and Sc is the collapse threshold. This equation can be used to 
derive the condition on the mean overdensity within a region of mass M in 
order for it to be ionized, 

Sm > Sb{M,z) = 5c- V2K{0[al,^ - a\M,z)]^'\ (170) 

where K{C,) — erfc~^(l — (^^^). Furlanetto et al. [143] showed how to construct 
the mass function of ionized regions from in analogy with the halo mass 
function (Press & Schechter 1974 [291]; Bond et al. 1991 [52]). The barrier in 
equation (170) is well approximated by a linear dependence on ct^, 

Sb ~ B{M) = Bo + Bia^iM), (171) 

in which case the mass function has an analytic solution (Sheth 1998 [332]), 



dluM 



5o 

■ exp 



cr(M) 



S2(M) 



2ct2(M) 



(172) 



where p is the mean mass density. This solution provides the comoving number 
density of ionized bubbles with mass in the range of (M, M + dM). The main 
difference of this result from the Press-Schechter mass function is that the 
barrier in this case becomes more difficult to cross on smaller scales because 
5b is a decreasing function of mass M. This gives bubbles a characteristic size. 
The size evolves with redshift in a way that depends only on Q and Mmin- 

One limitation of the above analytic model is that it ignores the non- 
local influence of sources on distant regions (such as voids) as well as the 
possible shadowing effect of intervening gas. Radiative transfer effects in the 
real Universe are inherently three-dimensional and cannot be fully captured 
by spherical averages as done in this model. Moreover, the value of Mmin is 
expected to increase in regions that were already ionized, complicating the ex- 
pectation of whether they will remain ionized later. The history of reionization 
could be complicated and non monotonic in individual regions, as described 
by Furlanetto & Loeb (2005) [144]. Finally, the above analytic formalism does 
not take the light propagation delay into account as we have done above in 
estimating the characteristic bubble size at the end of reionization. Hence this 
formalism describes the observed bubbles only as long as the characteristic 
bubble size is sufSciently small, so that the light propagation delay can be 
neglected compared to cosmic variance. The general effect of the light prop- 
agation delay on the power-spectrum of 21cm fluctuations was quantified by 
Barkana & Loeb (2005) [29]. 
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9.3 Separating the "Physics" from the "Astrophysics" of the 
Reionization Epoch with 21cm Fluctuations 

The 21cm signal can be seen from epochs during which the cosmic gas was 

largely neutral and deviated from thermal equilibrium with the cosmic mi- 
crowave background (CMB). The signal vanished at redshifts z > 200, when 
the residual fraction of free electrons after cosmological recombination kept 
the gas kinetic temperature, Tk, close to the CMB temperature, T^. But dur- 
ing 200 > z > 30 the gas cooled adiabatically and atomic collisions kept the 
spin temperature of the hyperfine level population below T^, so that the gas 
appeared in absorption [323, 226]. As the Hubble expansion continued to rar- 
efy the gas, radiative coupling of Tg to began to dominate and the 21cm 
signal faded. When the first galaxies formed, the UV photons they produced 
between the Lya and Lyman limit wavelengths propagated freely through 
the Universe, redshifted into the Lya resonance, and coupled Ts and Tk once 
again through the Wouthuysen-Field [388, 131] effect by which the two hy- 
perfine states are mixed through the absorption and re-emission of a Lya 
photon [237, 96]. Emission above the Lyman limit by the same galaxies ini- 
tiated the process of reionization by creating ionized bubbles in the neutral 
cosmic gas, while X-ray photons propagated farther and heated Tk above T^ 
throughout the Universe. Once Tg grew larger than T^, the gas appeared in 
21cm emission. The ionized bubbles imprinted a knee in the power spectrum 
of 21cm fluctuations [404], which traced the H I topology until the process of 
reionization was completed [143]. 

The various effects that determine the 21cm fluctuations can be separated 
into two classes. The density power spectrum probes basic cosmological pa- 
rameters and inflationary initial conditions, and can be calculated exactly in 
linear theory. However, the radiation from galaxies, both Lya radiation and 
ionizing photons, involves the complex, non- linear physics of galaxy forma- 
tion and star formation. If only the sum of all fluctuations could be measured, 
then it would be difficult to extract the separate sources, and in particular, 
the extraction of the power spectrum would be subject to systematic errors in- 
volving the properties of galaxies. Barkana & Loeb (2005) [28] showed that the 
unique three-dimensional properties of 21cm measurements permit a separa- 
tion of these distinct effects. Thus, 21cm fluctuations can probe astrophysical 
(radiative) sources associated with the first galaxies, while at the same time 
separately probing the physical (inflationary) initial conditions of the Uni- 
verse. In order to affect this separation most easily, it is necessary to measure 
the three-dimensional power spectrum of 21cm fluctuations. The discussion 
in this section follows Barkana & Loeb (2005) [28]. 
Spin temperature history 

As long as the spin-temperature T, is smaller than the CMB temperature 
T^ = 2.725(1 + z) K, hydrogen atoms absorb the CMB, whereas if Ts > T^ 
they emit excess flux. In general, the resonant 21cm interaction changes the 
brightness temperature of the CMB by [323, 237] U = t (T, - T^) /(I -t- z), 
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where the optical depth at a wavelength A = 21cm is 

ScX^hAionn 



327rA;Ts (1 + z) [dvr/dr. 



(173) 



where Hh is the number density of hydrogen, A\q = 2.85 x 10~^ s~ is the 

spontaneous emission coefficient, x^ii is the neutral hydrogen fraction, and 
dVr/dr is the gradient of the radial velocity along the line of sight with Vr being 
the physical radial velocity and r the comoving distance; on average dVr/dr = 
H{z)/{1 + z) where H is the Hubble parameter. The velocity gradient term 
arises because it dictates the path length over which a 21cm photon resonates 
with atoms before it is shifted out of resonance by the Doppler effect [341]. 

For the concordance set of cosmological parameters [348] , the mean bright- 
ness temperature on the sky at rcdshift z is 



Tfe = 28 mK 



T27 



(1 + ^) 
10 



1/2 



(174) 



where xhi is the mean neutral fraction of hydrogen. The spin temperature it- 
self is coupled to Tk through the spin-flip transition, which can be excited by 
collisions or by the absorption of Lya photons. As a result, the combination 
that appears in n becomes [131] (T, - T^)/r, = [a;tot/(l +a;tot)] (1 - T^/Tk), 
where x^ot = Xa + Xc is the sum of the radiative and coUisional thresh- 
old parameters. These parameters are Xa = 4:PaTi,/27AiQTry and Xc = 
4Ki_o(Tfe) nuTi^/ZAioT^, where Pa is the Lya scattering rate which is propor- 
tional to the Lya intensity, and ki_o is tabulated as a function of [11, 406]. 
The coupling of the spin temperature to the gas temperature becomes sub- 
stantial when Xtot > 1- 
Brightness temperature fluctuations 

Although the mean 21cm emission or absorption is difficult to measure due 
to bright foregrounds, the unique character of the ffuctuations in Tt, allows for 
a much easier extraction of the signal [154, 404, 259, 260, 314]. We adopt the 
notation 5a for the fractional fluctuation in quantity A (with a lone 5 denoting 
density perturbations). In general, the fluctuations in Tjj can be sourced by 
fluctuations in gas density (5), Lya flux (through 5x^) neutral fraction (5xhi), 
radial velocity gradient {5d^y^), and temperature, so we flnd 



1 + ^ 

Xtot 
+ (7a - 1) 



Xtot 



+ - 



Xc dlog(Ki_o) 



Tk-T^ Xtot rflog(Tfc) 



(175) 



where the adiabatic index is 7a = 1 -|- {St^/S), and we deflne Xtot = (1 + 
a;tot)a;tot- Taking the Fourier transform, we obtain the power spectrum of 
each quantity; e.g., the total power spectrum Pj-^ is defined by 



{5tM5tM) = (27r)^5^(ki +k2)PT,(ki) 



(176) 
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where Sr^i^) is the Fourier transform of Sx,,, k is the comoving wavevector, 
6^ is the Dirac delta function, and (• • • ) denotes an ensemble average. In this 
analysis, we consider scales much bigger than the characteristic bubble size 
and the early phase of reionization (when Sx^i « 1); so that the fluctuations 
Sx-iij are also much smaller than unity. For a more general treatment, see 
McQuinn et al. (2005) [250]. 

The separation of powers 

The fluctuation St^ consists of a number of isotropic sources of fluctua- 
tions plus the peculiar velocity term —6^^^^. Its Fourier transform is simply 
proportional to that of the density field [191, 41], 

~Sd^v^ = -fJ-^l (177) 

where /i = cos 0k in terms of the angle 9k of k with respect to the line of 
sight. The /x^ dependence in this equation results from taking the radial (i.e., 
line-of-sight) component (oc /x) of the peculiar velocity, and then the radial 
component (oc /i) of its gradient. Intuitively, a high-density region possesses a 
velocity infall towards the density peak, implying that a photon must travel 
further from the peak in order to reach a fixed relative redshift, compared with 
the case of pure Hubble expansion. Thus the optical depth is always increased 
by this effect in regions with 5 > 0. This phenomenon is most properly termed 
velocity compression. 

We therefore write the fluctuation in Fourier space as 

5T,(k) = fi^S{k) + f3S{k) + (5,ad(k) , (178) 

where we have defined a coefficient /3 by collecting all terms oc 5 in Eq. (175), 
and have also combined the terms that depend on the radiation fields of 
Lya photons and ionizing photons, respectively. We assume that these ra- 
diation fields produce isotropic power spectra, since the physical processes 
that determine them have no preferred direction in space. The total power 
spectrum is 

PT,(k) = fx^Psik) + 2ti^[(3Ps{k) + Ps.rM] + 

[P'^Psik) + P.ad(fc) + 2pPs.reA{k)] , (179) 

where we have defined the power spectrum P^.rad as the Fourier transform of 
the cross-correlation function, 

6.rad(r) = (5(ri) 5,ad(ri + r)) . (180) 

We note that a similar anisotropy in the power spectrum has been pre- 
viously derived in a different context, i.e., where the use of galaxy redshifts 
to estimate distances changes the apparent line-of-sight density of galaxies in 
redshift surveys [191, 219, 178, 133]. However, galaxies arc intrinsically com- 
plex tracers of the underlying density field, and in that case there is no analog 
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to the method that we demonstrate below for separating in 21cm fluctuations 
the effect of initial conditions from that of later astrophysical processes. 

The velocity gradient term has also been examined for its global eS'ect on 
the sky-averaged power and on radio visibilities [366, 41]. The other sources 
of 21cm perturbations arc isotropic and would produce a power spectrum 
Pt^ (fc) that could be measured by averaging the power over spherical shells 
in k space. In the simple case where /3 = 1 and only the density and velocity 
terms contribute, the velocity term increases the total power by a factor of 
((1 + n'^Y) = 1.87 in the spherical average. However, instead of averaging 
the signal, we can use the angular structure of the power spectrum to greatly 
increase the discriminatory power of 21cm observations. We may break up each 
spherical shell in k space into rings of constant /z and construct the observed 

[k, n). Considering Eq. (179) as a polynomial in fi, i.e., ii^Pni+fi^P^-i +PijO , 
wc see that the power at just three values of /x is required in order to separate 
out the coefficients of 1, /x^, and /x^ for each k. 

If the velocity compression were not present, then only the /i-independent 
term (times T^) would have been observed, and its separation into the five 
components {Tj,, jS, and three power spectra) would have been difficult and 
subject to degeneracies. Once the power has been separated into three parts, 
however, the /i^ coefficient can be used to measure the density power spec- 
trum directly, with no interference from any other source of fluctuations. 
Since the overall amplitude of the power spectrum, and its scaling with 
redshift, arc well determined from the combination of the CMB tempera- 
ture fluctuations and galaxy surveys, the amplitude of P^4 directly deter- 
mines the mean brightness temperature Tb on the sky, which measures a 
combination of and Xfji at the observed redshift. McQuinn et al. (2005) 
[250] analysed in detail the parameters that can be constrained by upcom- 
ing 21cm experiments in concert with future CMB experiments such as 
Planck (http://www.rssd.esa.int/indcx.php?projcct=PLANCK). Once Ps{k) 
has been determined, the coefficients of the /x^ term and the /x-independent 
term must be used to determine the remaining unknowns, /3, Ps-ra,d{k), and 
Piad{k). Since the coefficient {3 is independent of fc, determining it and thus 
breaking the last remaining degeneracy requires only a weak additional as- 
sumption on the behavior of the power spectra, such as their asymptotic 
behavior at large or small scales. If the measurements cover Nk values of 
wavenumber k, then one wishes to determine 2Nk + 1 quantities based on 2Nh 
measurements, which should not cause significant degeneracies when ^ 1. 
Even without knowing /3, one can probe whether some sources of Prad(fc) 
are uncorrclatcd with d; the quantity Pun-s{k) = P^o — P^2/(4P^4) equals 
-Prad — Ps rad/-Ps, which receives no contribution from any source that is a 
linear functional of the density distribution (see the next subsection for an 
example) . 
Specific epoclis 

At z ^ 35, collisions are effective due to the high gas density, so one can 
measure the density power spectrum [226] and the redshift evolution of nm, 
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T^, and T^. At z < 35, collisions become ineffective but the first stars pro- 
duce a cosmic background of Lya photons (i.e. photons that redshift into the 
Lyof resonance) that couples Tg to T^. During the period of initial Lya cou- 
pling, fluctuations in the Lya flux translate into fluctuations in the 21cm 
brightness [30]. This signal can be observed from z ~ 25 until the Lya cou- 
phng is completed (i.e., Xtot ^ 1) at z ~ 15. At a given redshift, each atom 
sees Lya photons that were originally emitted at earlier times at rest-frame 
wavelengths between Lya and the Lyman limit. Distant sources are time re- 
tarded, and since there are fewer galaxies in the distant, earlier Universe, each 
atom sees sources only out to an apparent source horizon of 100 comoving 
Mpc at z ^ 20. A significant portion of the flux comes from nearby sources, 
because of the 1/r^ decline of flux with distance, and since higher Lyman 
series photons, which are degraded to Lya photons through scattering, can 
only be seen from a small redshift interval that corresponds to the wavelength 
interval between two consecutive atomic levels. 




Wavenumber [1/Mpc] 

Fig. 63. Observable power spectra during the period of initial Lya coupling. Upper 
panel: Assumes adiabatic cooling. Lower panel: Assumes pre-heating to 500 K by 
X-ray sources. Shown are P^i = Ps (sohd curves), P^2 (short-dashed curves), and 
Pun-s (long-dashed curves), as well as for comparison 2/3Ps (dotted curves). 

There are two separate sources of fluctuations in the Lya flux [30]. The 
first is density inhomogeneities. Since gravitational instability proceeds faster 
in overdense regions, the biased distribution of rare galactic halos fluctuates 
much more than the global dark matter density. When the number of sources 
seen by each atom is relatively small, Poisson fluctuations provide a second 
source of fluctuations. Unlike typical Poisson noise, these fluctuations are cor- 
related between gas elements at different places, since two nearby elements see 
many of the same sources. Assuming a scale-invariant spectrum of primordial 
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density fluctuations, and that = 1 is produced at 2; = 20 by galaxies in 
dark matter halos where the gas cools cfScicntly via atomic cooling, Figure 
63 shows the predicted observable power spectra. The figure suggests that 
/? can be measured from the ratio P^-i/P^i at fc > 1 Mpc~^, allowing the 
density-induced fluctuations in flux to be extracted from P^,2. while only the 
Poisson fluctuations contribute to Pun-5- Each of these components probes 
the number density of galaxies through its magnitude, and the distribution 
of source distances through its shape. Measurements at fc > 100 Mpc"-'^ can 
independently probe because of the smoothing effects of the gas pressure 
and the thermal width of the 21cm line. 

After Lya coiipling and X-ray heating are both completed, reionization 
continues. Since /3 = 1 and ^ T^, the normalization of directly 
measures the mean neutral hydrogen fraction, and one can separately probe 
the density fluctuations, the neutral hydrogen fluctuations, and their cross- 
correlation. 

Fluctuations on large angular scales 

Full-sky observations must normally be analyzed with an angular and ra- 
dial transform [143, 314, 41], rather than a Fourier transform which is simpler 
and yields more directly the underlying 3D power spectrum [259, 260]. The 
21cm brightness fluctuations at a given redshift - corresponding to a comov- 
ing distance ro from the observer - can be expanded in spherical harmonics 
with expansion coefficients aim{i^), where the angular power spectrum is 



with Gi{x) = Ji{x) + (/? — l)ji{x) and J; (a;) being a linear combination of 
spherical Bessel functions [41]. 

In an angular transform on the sky, an angle of 9 radians translates to a 
spherical multipolc I ^ 3.5/6. For measurements on a screen at a comoving 
distance ro, a multipole I normally measures 3D power on a scale oi ~ 
9ro ~ 35/Z Gpc for I » 1, since ro ~ 10 Gpc at z > 10. This estimate fails 
at I < 100, however, when we consider the sources of 21cm fluctuations. The 
angular projection implied in C; involves a weighted average (Eq. 181) that 
favors large scales when I is small, but density fluctuations possess little large- 
scale power, and the C/ are dominated by power around the peak of kPs{k), 
at a few tens of comoving Mpc. 

Figure 64 shows that for density and velocity fluctuations, even the I = 1 
multipole is affected by power at fc"-' > 200 Mpc only at the 2% level. Due to 
the small number of large angular modes available on the sky, the expectation 
value of Ci cannot be measured precisely at small I. Figure 64 shows that 
this precludes new information from being obtained on scales > 130 Mpc 
using angular structure at any given redshift. Fluctuations on such scales may 
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Fig. 64. Effect of large-scale power on the angular power spectrum of 21cm 
anisotropies on the sky. This example shows the power from density fluctuations 
and velocity compression, assuming a warm IGM at z = 12 with Ts = S> T-y. 
Shown is the % change in Ci if we were to cut off the power spectrum above l/fc of 
200, 180, 160, 140, 120, and 100 Mpc (top to bottom). Also shown for comparison 
is the cosmic variance for averaging in bands of ZiZ ~ i (dashed lines). 

be measurable using a range of redshifts, but the required Az > 1 at z ~ 10 
implies significant difficulties with foreground subtraction and with the need 
to account for time evolution. 

10 Major Challenge for Future Theoretical Research: 

radiative transfer during reionization requires a large dynamic 
range, challenging the capabilities of existing simulation codes 

Observations of the cosmic microwave background [348] have confirmed the 
notion that the present large-scale structure in the Universe originated from 
small-amplitude density fluctuations at early cosmic times. Due to the natural 
instability of gravity, regions that were denser than average collapsed and 
formed bound halos, first on small spatial scales and later on larger and larger 
scales. At each snapshot of this cosmic evolution, the abundance of collapsed 
halos, whose masses are dominated by cold dark matter, can be computed 
from the initial conditions using numerical simulations and can be understood 
using approximate analytic models [292, 52]. The common understanding of 
galaxy formation is based on the notion that the constituent stars formed out 
of the gas that cooled and subsequently condensed to high densities in the 
cores of some of these halos [379]. 

The standard analytic model for the abundance of halos [292, 52] considers 
the small density fluctuations at some early, initial time, and attempts to 



First Light 137 



predict the number of halos that will form at some later time corresponding 

to a rcdshift z. First, the fluctuations arc extrapolated to the present time 
using the growth rate of linear fluctuations, and then the average density 
is computed in spheres of various sizes. Whenever the overdensity (i.e., the 
density perturbation in units of the cosmic mean density) in a sphere rises 
above a critical threshold (5c (2), the corresponding region is assumed to have 
collapsed by redshift 0, forming a halo out of all the mass that had been 
included in the initial spherical region. In analyzing the statistics of such 
regions, the model separates the contribution of large-scale modes from that 
of small-scale density fluctuations. It predicts that galactic halos will form 
earlier in regions that arc overdense on large scales [190, 19, 97, 258], since 
these regions already start out from an enhanced level of density, and small- 
scale modes need only supply the remaining perturbation necessary to reach 
(5c (z). On the other hand, large-scale voids should contain a reduced number of 
halos at high redshift. In this way, the analytic model describes the clustering 
of massive halos. 

As gas falls into a dark matter halo, it can fragment into stars only if its 
virial temperature is above 10**K for cooling mediated by atomic transitions 
[or ~ 500 K for molecular H2 cooling; see Fig. 20]. The abundance of dark 
matter halos with a virial temperature above this cooling threshold declines 
sharply with increasing redshift due to the exponential cutoff in the abun- 
dance of massive halos at early cosmic times. Consequently, a small change 
in the collapse threshold of these rare halos, due to mild inhomogeneities on 
much larger spatial scales, can change the abundance of such halos dramati- 
cally. Barkana & Loeb (2004) [27] have shown that the modulation of galaxy 
formation by long wavelength modes of density fluctuations is therefore am- 
plified considerably at high redshift; the discussion in this section follows their 
analysis. 

Ampliflcation of Density Fluctuations 

Galaxies at high redshift are believed to form in all halos above some 
minimum mass Mmin that depends on the eflaciency of atomic and molecular 
transitions that cool the gas within each halo. This makes useful the standard 
quantity of the collapse fraction -Fcoi(Mmin), which is the fraction of mass in 
a given volume that is contained in halos of individual mass Mmin or greater 
(see Fig. 13). If wc set Afmin to be the minimum halo mass in which efficient 
cooling processes are triggered, then -Fed (Mmin) is the fraction of all baryons 
that reside in galaxies. In a large-scale region of comoving radius R with a 
mean overdensity bji, the standard result is [52] 



-Fcoi(Mmin) = erfc 



^c(^) - 



V2 [5(i?min) - S(R)\ 



(182) 



where S(R) = cr'^iR) is the variance of fluctuations in spheres of radius R, and 

5'(i?inin) is the variance in spheres of radius -Rmin corresponding to the region 
at the initial time that contained a mass Mmin- In particular, the cosmic mean 
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value of the collapse fraction is obtained in the limit of i? ^ c» by setting 

6fj and S{R) to zero in this expression. Throughout this section we shall 
adopt this standard model, known as the extended Press-Schechter model. 
Whenever we consider a cubic region, we will estimate its halo abundance 
by applying the model to a spherical region of equal volume. Note also that 
we will consistently quote values of comoving distance, which equals physical 
distance times a factor of {1 + z). 

At high redshift, galactic halos are rare and correspond to high peaks in 
the Gaussian probability distribution of initial fluctuations. A modest change 
in the overall density of a large region modulates the threshold for high peaks 
in the Gaussian density field, so that the number of galaxies is exponentially 
sensitive to this modulation. This amplification of large-scale modes is respon- 
sible for the large statistical fluctuations that we find. 

In numerical simulations, periodic boundary conditions arc usually as- 
sumed, and this forces the mean density of the box to equal the cosmic mean 
density. The abundance of halos as a function of mass is then biased in such a 
box (see Fig. 65) , since a similar region in the real Universe will have a distri- 
bution of different overdensities Sr. At high redshift, when galaxies correspond 
to high peaks, they are mostly found in regions with an enhanced large-scale 
density. In a periodic box, therefore, the total number of galaxies is artificially 
reduced, and the relative abundance of galactic halos with different masses is 
artificially tilted in favor of lower-mass halos. Let us illustrate these results 
for two sets of parameters, one corresponding to the first galaxies and early 
reionization [z = 20) and the other to the current horizon in observations 
of galaxies and late reionization {z = 7). Let us consider a resolution equal 
to that of state-of-the-art cosmological simulations that include gravity and 
gas hydrodynamics. Specifically, let us assume that the total number of dark 
matter particles in the simulation is N = 324^, and that the smallest halo 
that can form a galaxy must be resolved into 500 particles; [349] showed that 
this resolution is necessary in order to determine the star formation rate in an 
individual halo reliably to within a factor of two. Therefore, if we assume that 
halos that cool via molecular hydrogen must be resolved at z = 20 (so that 
Afmin = 7 X IO^Mq), and only those that cool via atomic transitions must 
be resolved at ^ = 7 (so that Mmin = W^Mq), then the maximum box sizes 
that can currently be simulated in hydrodynamic comological simulations are 
^box = 1 Mpc and /box = 6 Mpc at these two redshifts, respectively. 

At each redshift we only consider cubic boxes large enough so that the 
probability of forming a halo on the scale of the entire box is negligible. In 
this case, Sr is Gaussian distributed with zero mean and variance S{R), since 
the no- halo condition ^/ S{R) Sc{z) implies that at redshift z the per- 
turbation on the scale R is still in the linear regime. We can then calculate 
the probability distribution of collapse fractions in a box of a given size (see 
Fig.66). This distribution corresponds to a real variation in the fraction of 
gas in galaxies within different regions of the Universe at a given time. In a 
numerical simulation, the assumption of periodic boundary conditions elimi- 
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Fig. 65. Bias in the halo mass distribution in simulations. Shown is the amount 
of mass contained in all halos of individual mass Mmin or greater, expressed as a 
fraction of the total mass in a given volume. This cumulative fraction Fcoi(Mi„in) is 
illustrated as a function of the minimum halo mass Mmin. We consider two cases of 
redshift and simulation box size, namely z — 7, /box = 6 Mpc (upper curves), and 
z = 20, Zbox = 1 Mpc (lower curves). At each redshift, we compare the true average 
distribution in the Universe (dotted curve) to the biased distribution (solid curve) 
that would be measured in a simulation box with periodic boundary conditions (for 
which Zr is artificially set to zero). 



nates the large-scale modes that cause this cosmic scatter. Note that Poisson 
fluctuations in the number of halos within the box would only add to the 
scatter, although the variations we have calculated are typically the domi- 
nant factor. For instance, in our two standard examples, the mean expected 
number of halos in the box is 3 at 2; = 20 and 900 at 2; = 7, resulting in 
Poisson fluctuations of a factor of about 2 and 1.03, respectively, compared 
to the clustering-induced scatter of a factor of about 16 and 2 in these two 
cases. 

Within the extended Press-Schechter model, both the numerical bias and 
the cosmic scatter can be simply described in terms of a shift in the redshift 
(see Fig. 67). In general, a region of radius R with a mean overdensity 5 a will 
contain a different collapse fraction than the cosmic mean value at a given 
redshift z. However, at some wrong redshift z -f Az this small region will 
contain the cosmic mean collapse fraction at z. At high redshifts (z > 3), this 
shift in redshift was derived by Barkana & Loeb [27] from equation (182) [and 
was already mentioned in Eq. (168)] 



Ziz = ^-(1 + z) X 



1- Wl - 



S{R) 

S{Rmin) 



(183) 
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Fig. 66. Probability distribution within a small volume of the total mass frac- 
tion in galactic halos. The normalized distribution of the logarithm of this fraction 
Jcoi(Afmin) is showu for two cases: z — 7, Zbox = 6 Mpc, Mmin = 10* Mq (upper 
panel), and z = 20, Zbox = 1 Mpc, Mmin = 7 x IQ^Mq (bottom panel). In each case, 
the value in a periodic box {5r = 0) is shown along with the value that would be 
expected given a plus or minus 1 — a fluctuation in the mean density of the box 
(dashed vertical lines). Also shown in each case is the mean value of -F'coi(Mmin) 
averaged over large cosmological volumes (solid vertical line). 



where 5q = 5c{z)/{l + z) is approximately constant at high redshifts [283], 
and equals 1.28 for the standard cosmological parameters (with its deviation 
from the Einstein-de Sitter value of 1.69 resulting from the existence of a 
cosmological constant). Thus, in our two examples, the bias is -2.6 at 2 = 20 
and -0.4 at z = 7, and the one-sided 1 — cr scatter is 2.4 at z = 20 and 1.2 at 
z = 7. 

Matching Numerical Simulations 

Next we may develop an improved model that fits the results of numerical 
simulations more accurately. The model constructs the halo mass distribution 
(or mass function); cumulative quantities such as the collapse fraction or the 
total number of galaxies can then be determined from it via integration. We 
first define f{Sc{z), S) dS to be the mass fraction contained at z within halos 
with mass in the range corresponding to S to S + dS. As derived earlier, the 
Press-Schechter halo abundance is 



dn po 
dM ^ M 



dS 



dM 



fiSc{z),S), (184) 



First Light 141 



1 E — I I I 




JQ-3 I I I I I I I I I I I I 

0.1 1 10' 10= 10= 

Box size [Mpc] 

Fig. 67. Cosmic scatter and numerical bias, expressed as the change in redshift 
needed to get the correct cosmic mean of the collapse fraction. The plot shows the 
1 — a scatter (about the biased value) in the redshift of reionization, or any other 
phenomenon that depends on the mass fraction in galaxies (bottom panel), as well 
as the redshift bias [expressed as a fraction of (1 + z)] in periodic simulation boxes 
(upper panel). The bias is shown for Mmin = 7 x 10''' Af© (solid curve), Afmin ~ 
10* M© (dashed curve), and Mmin = 3 x 10^^ M© (dotted curve). The bias is always 
negative, and the plot gives its absolute value. When expressed as a shift in redshift, 
the scatter is independent of Mmin. 



where dn is the comoving number density of halos with masses in the range 
M to M + dM, and 

(185) 

where v = 5c{z)/y/S is the number of standard deviations that the critical 
collapse overdensity represents on the mass scale M corresponding to the 
variance S. 

However, the Press-Schechter mass function fits numerical simulations only 
roughly, and in particular it substantially underestimates the abundance of 
the rare halos that host galaxies at high redshift. The halo mass function of 
[333] [see also [334]] adds two free parameters that allow it to fit numerical 
simulations much more accurately [188]. These N-body simulations followed 
very large volumes at low redshift, so that cosmic scatter did not compromise 
their accuracy. The matching mass function is given by 



exp 
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fsTi6c{z),S) 



1 + 





■ aV2- 


exp 


2 



(186) 



(aV2)9' 

0.3, and where normalization 



with best-fit parameters [335] a' = 0.75 and q' 

to unity is ensured by taking A' = 0.322. 

In order to calculate cosmic scatter we must determine the biased halo 
mass function in a given volume at a given mean density. Within the extended 
Press-Schechter model [52], the halo mass distribution in a region of comoving 
radius R with a mean ovcrdcnsity Sn is given by 



/bias-ps(<5c(2), Sr, R, S) = fps{6c{z) -5r,S- S{R)) . 



(187) 



The corresponding collapse fraction in this case is given simply by eq. (182). 
Despite the relatively low accuracy of the Press-Schechter mass function, the 
relative change is predicted rather accurately by the extended Press-Schechter 
model. In other words, the prediction for the halo mass function in a given 
volume compared to the cosmic mean mass function provides a good fit to 
numerical simulations over a wide range of parameters [258, 77] . 

For the improved model (derived in [27]), we adopt a hybrid approach that 
combines various previous models with each applied where it has been found 
to closely match numerical simulations. We obtain the halo mass function 
within a restricted volume by starting with the Sheth-Torme formula for the 
cosmic mean mass function, and then adjusting it with a relative correction 
based on the extended Press-Schechter model. In other words, we set 



fbi^i6ciz),dR,R, S) = 



fps{6ciz),S) 



(188) 



As noted, this model is based on fits to simulations at low redshifts, but we can 
check it at high redshifts as well. Figure 68 shows the number of galactic halos 
at 2; ~ 15 — 30 in two numerical simulations run by [402], and our predictions 
given the cosmological input parameters assumed by each simulation. The 
close fit to the simulated data (with no additional free parameters) suggests 
that our hybrid model (solid linos) improves on the extended Press-Schechter 
model (dashed lines), and can bo used to calculate accurately the cosmic 
scatter in the number of galaxies at both high and low redshifts. The simulated 
data significantly deviate from the expected cosmic mean [eq. (186), shown by 
the dotted line] , due to the artificial suppression of large-scale modes outside 
the simulated box. 

As an additional example, we consider the highest-resolution first star 
simulation [5], which used /box = 128 kpc and A'/,nin = 7 x 10^ Mq. The first 
star forms within the simulated volume when the first halo of mass Mmin or 
larger collapses within the box. To compare with the simulation, we predict 
the redshift at which the probability of finding at least one halo within the box 
equals 50%, accounting for Poisson fiuctuations. We find that if the simulation 
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Redshift 

Fig. 68. Halo mass lunction at high redshift in a 1 Mpc box at the cosmic mean 
density. The prediction (solid lines) of the hybrid model of Barkana & Loeb (2004) 
[27] is compared with the number of halos above mass 7 x 10^ M0 measured in the 
simulations of [402] [data points are taken from their Figure 5] . The cosmic mean of 
the halo mass function (dotted lines) deviates significantly from the simulated values, 
since the periodic boundary conditions within the finite simulation box artificially 
set the amplitude of large-scale modes to zero. The hybrid model starts with the 
Sheth-Tormen mass function and applies a correction based on the extended Press- 
Schechter model; in doing so, it provides a better fit to numerical simulations than 
the pure extended Press-Schechter model (dashed lines) used in the previous figures. 
We consider two sets of cosmological parameters, the scale-invariant ylCDM model 
of [402] (upper curves), and their running scalar index (RSI) model (lower curves). 

formed a population of halos corresponding to the correct cosmic average [as 
given by eq. (186)], then the first star should have formed already dX z — 24.0. 
The first star actually formed in the simulation box only at z — 18.2 [5]. Using 
eq. (188) we can account for the loss of large-scale modes beyond the periodic 
box, and predict a first star at z = 17.8, a close match given the large Poisson 
fluctuations introduced by considering a single galaxy within the box. 

The artificial bias in periodic simulation boxes can also be seen in the 
results of extensive numerical convergence tests carried out by [349]. They 
presented a large array of numerical simulations of galaxy formation run in 
periodic boxes over a wide range of box size, mass resolution, and redshift. In 
particular, we can identify several pairs of simulations where the simulations 
in each pair have the same mass resolution but different box sizes; this allows 
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us to separate the effect of large-scale numerical bias from the effect of having 
poorly-resolved individual halos. 

Implications 

(i) The nature of reionization 

A variety of papers in the literature [13, 138, 330, 171, 152, 23, 224] main- 
tain that reionization ended with a fast, simultaneous, overlap stage through- 
out the Universe. This view has been based on simple arguments and has been 
supported by numerical simulations with small box sizes. The underlying idea 
was that the ionized hydrogen (H II ) regions of individual sources began to 
overlap when the typical size of each H II bubble became comparable to the 
distance between nearby sources. Since these two length scales were compa- 
rable at the critical moment, there is only a single timescale in the problem 
- given by the growth rate of each bubble - and it determines the transition 
time between the initial overlap of two or three nearby bubbles, to the final 
stage where dozens or hundreds of individual sources overlap and produce 
large ionized regions. Whenever two ionized bubbles were joined, each point 
inside their common boundary became exposed to ionizing photons from both 
sources, reducing the neutral hydrogen fraction and allowing ionizing photons 
to travel farther before being absorbed. Thus, the ionizing intensity inside 
H II regions rose rapidly, allowing those regions to expand into high-density 
gas that had previously recombined fast enough to remain neutral when the 
ionizing intensity had been low. Since each bubble coalescence accelerates the 
process, it has been thought that the overlap phase has the character of a 
phase transition and occurs rapidly. Indeed, the simulations of reionization 
[152] found that the average mean free path of ionizing photons in the simu- 
lated volume rises by an order of magnitude over a redshift interval Az = 0.05 
at 2; = 7. 

These results imply that overlap is still expected to occur rapidly, but only 
in localized high-density regions, where the ionizing intensity and the mean 
free path rise rapidly even while other distant regions are still mostly neutral. 
In other words, the size of the bubble of an individual source is about the same 
in different regions (since most halos have masses just above Mmin), but the 
typical distance between nearby sources varies widely across the Universe. 
The strong clustering of ionizing sources on length scales as large as 30- 
100 Mpc introduces long timescales into the reionization phase transition. 
The sharpness of overlap is determined not by the growth rate of bubbles 
around individual sources, but by the ability of large groups of sources within 
overdense regions to deliver ionizing photons into large underdense regions. 

Note that the recombination rate is higher in overdense regions because of 
their higher gas density. These regions still reionize first, though, despite the 
need to overcome the higher recombination rate, since the number of ionizing 
sources in these regions is increased even more strongly as a result of the 
dramatic amplification of large-scale modes discussed earlier. 

(ii) Limitations of current simulations 



First Light 145 



The shortcomings of current simulations do not amount simply to a shift 
of ^ 10% in redshift and the elimination of scatter. The effect mentioned 
above can be expressed in terms of a shift in redshift only within the context 
of the extended Press-Schechter model, and only if the total mass fraction in 
galaxies is considered and not its distribution as a function of galaxy mass. 
The halo mass distribution should still have the wrong shape, resulting from 
the fact that Az depends on Mmin- A self-contained numerical simulation 
must directly evolve a very large volume. 

Another reason that current simulations are limited is that at high redshift, 
when galaxies are still rare, the abundance of galaxies grows rapidly towards 
lower redshift. Therefore, a ^ 10% relative error in redshift implies that at 
any given redshift around z ~ 10-20, the simulation predicts a halo mass 
function that can be off by an order of magnitude for halos that host galaxies 
(see Fig. 68). This large underestimate suggests that the first generation of 
galaxies formed significantly earlier than indicated by recent simulations. An- 
other element missed by simulations is the large cosmic scatter. This scatter 
can fundamentally change the character of any observable process or feedback 
mechanism that depends on a radiation background. Simulations in periodic 
boxes eliminate any large-scale scatter by assuming that the simulated volume 
is surrounded by identical periodic copies of itself. In the case of rcionization, 
for instance, current simulations neglect the collective effects described above, 
whereby groups of sources in overdense regions may influence large surround- 
ing underdense regions. In the case of the formation of the first stars due to 
molecular hydrogen cooling, the effect of the soft ultraviolet radiation from 
these stars, which tends to dissociate the molecular hydrogen around them 
[170, 303, 272], must be reassessed with cosmic scatter included. 
(in) Observational consequences 

The spatial fluctuations that we have calculated also affect current and 
future observations that probe rcionization or the galaxy population at high 
redshift. For example, there are a large number of programs searching for 
galaxies at the highest accessible redshifts (6.5 and beyond) using their strong 
Lya emission [184, 301, 244, 202]. These programs have previously been justi- 
fied as a search for the rcionization redshift, since the intrinsic emission should 
be absorbed more strongly by the surrounding IGM if this medium is neutral. 
For any particular source, it will be hard to clearly recognize this enhanced 
absorption because of uncertainties regarding the properties of the source and 
its radiative and gravitational effects on its surroundings [24, 26, 312]. How- 
ever, if the luminosity function of galaxies that emit Lya can be observed, 
then faint sources, which do not significantly affect their environment, should 
be very strongly absorbed in the era before rcionization. Rcionization can 
then be detected statistically through the sudden jump in the number of faint 
sources [246]. The above results alter the expectation for such observations. In- 
deed, no sharp "rcionization redshift" is expected. Instead, a Lya luminosity 
function assembled from a large area of the sky will average over the cosmic 
scatter of Az ~ 1-2 between different regions, resulting in a smooth evolu- 
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tion of the luminosity function over this redshift range. In addition, such a 

survey may be biased to give a relatively high redshift, since only the most 
massive galaxies can be detected, and as we have shown, these galaxies will be 
concentrated in overdense regions that will also get reionized relatively early. 

The distribution of ionized patches during reionization will likely be probed 
by future observations, including small-scale anisotropics of the cosmic mi- 
crowave background photons that are rescattered by the ionized patches 
[8, 162, 313], and observations of 21 cm emission by the spin-flip transition of 
the hydrogen in neutral regions [366, 75, 142]. Previous analytical and numer- 
ical estimates of these signals have not included the collective effects discussed 
above, in which rare groups of massive galaxies may reionize large surround- 
ing areas. The transfer of photons across large scales will likely smooth out 
the signal even on scales significantly larger than the typical size of an H II 
bubble due to an individual galaxy. Therefore, even the characteristic angular 
scales that are expected to show correlations in such observations must be 
reassessed. 

The cosmic scatter also affects observations in the present-day Universe 
that depend on the history of reionization. For instance, photoionization heat- 
ing suppresses the formation of dwarf galaxies after reionization, suggesting 
that the smallest galaxies seen today may have formed prior to reionization 
[73, 344, 37]. Under the popular view that assumed a sharp end to reioniza- 
tion, it was expected that denser regions would have formed more galaxies by 
the time of reionization, possibly explaining the larger relative abundance of 
dwarf galaxies observed in galaxy clusters compared to lower-density regions 
such as the Local Group of galaxies [369, 38]. The above results undercut 
the basic assumption of this argument and suggest a different explanation. 
Reionization occurs roughly when the number of ionizing photons produced 
starts to exceed the number of hydrogen atoms in the surrounding IGM. If 
the processes of star formation and the production of ionizing photons are 
equally efficient within galaxies that lie in different regions, then reionization 
in each region will occur when the collapse fraction reaches the same critical 
value, even though this will occur at different times in different regions. Since 
the galaxies responsible for reionization have the same masses as present-day 
dwarf galaxies, this estimate argues for a roughly equal abundance of dwarf 
galaxies in all environments today. This simple picture is, however, modifled 
by several additional effects. First, the recombination rate is higher in over- 
dense regions at any given time, as discussed above. Furthermore, reionization 
in such regions is accomplished at an earlier time when the recombination 
rate was higher even at the mean cosmic density; therefore, more ionizing 
photons must be produced in order to compensate for the enhanced recom- 
bination rate. These two effects combine to make overdense regions reionize 
at a higher value of F^oi than underdense regions. In addition, the overdense 
regions, which reionize first, subsequently send their extra ionizing photons 
into the surrounding underdense regions, causing the latter to reionize at an 
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even lower F^oi- Thus, a higher abundance of dwarf galaxies today is indeed 
expected in the ovcrdcnsc regions. 

The same basic effect may be even more critical for understanding the 
properties of large-scale voids, 10-30 Mpc regions in the present-day Universe 
with an average mass density that is well below the cosmic mean. In order to 
predict their properties, the first step is to consider the abundance of dark mat- 
ter halos within them. Numerical simulations show that voids contain a lower 
relative abundance of rare halos [249, 82, 39], as expected from the raising of 
the collapse threshold for halos within a void. On the other hand, simulations 
show that voids actually place a larger fraction of their dark matter content 
in dwarf halos of mass below lG^"A'/0 [157]. This can be understood within 
the extended Press- Schechter model. At the present time, a typical region in 
the Universe fills halos of mass lO^^M© and higher with most of the dark 
matter, and very little is left over for isolated dwarf halos. Although a large 
number of dwarf halos may have formed at early times in such a region, the 
vast majority later merged with other halos, and by the present time they 
survive only as substructure inside much larger halos. In a void, on the other 
hand, large halos are rare even today, implying that most of the dwarf halos 
that formed early within a void can remain as isolated dwarf halos till the 
present. Thus, most isolated dwarf dark matter halos in the present Universe 
should be found within large-scale voids [25] . 

However, voids are observed to be rather deficient in dwarf galaxies as 
well as in larger galaxies on the scale of the Milky Way mass of ^ 
[200, 120, 285]. A deficit of large galaxies is naturally expected, since the total 
mass density in the void is unusually low, and the fraction of this already 
low density that assembles in large halos is further reduced relative to higher- 
density regions. The absence of dwarf galaxies is harder to understand, given 
the higher relative abundance expected for their host dark matter halos. The 
standard model for galaxy formation may be consistent with the observations 
if some of the dwarf halos are dark and do not host stars. Large numbers of 
dark dwarf halos may be produced by the effect of reionization in suppressing 
the infall of gas into these halos. Indeed, exactly the same factors considered 
above, in the discussion of dwarf galaxies in clusters compared to those in small 
groups, apply also to voids. Thus, the voids should reionize last, but since they 
are most strongly affected by ionizing photons from their surroundings (which 
have a higher density than the voids themselves), the voids should reionize 
when the abundance of galaxies within them is relatively low. 
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